SR_SB_Veg_Analysis

About

Seed rain, seed bank, and vegetational dynamics of north-central Missouri remnant and reconstructed tallgrass prairies

Background:

Author: Katherine Wynne (Wynnekat@msu.edu)

Co-authors: E. Eyerly, L. Sullivan

Created: 14 June 2023

Files:

Cleaned_Spring_Cover_Data_July_2019.xlsx - contains spring p/a cover data
Cleaned_Fall_Cover_Data_September_2019.xlsx - contains fall % cover data
Seed_Rain_Data_2019_FINAL.xlsx - contains seed rain data
Seed_Bank_2020_Cleaned.xlsx - contains seed bank data
SR_Traits_Dataset.xlsx - contains trait data (work in progress)

R Version: R 4.3.2 (2023-10-31) – “Eye Holes”

RStudio Version: RStudio 2023.06.1+524 “Mountain Hydrangea”

Package Version:

Last updated: 23 February 2024

—–

Setup

Libraries

### --- Dataset importing and exporting -----

# Import in excel files
library(readxl)

# Export dataframes to excel files

library(writexl)

### ---- Data manipulation and visualization ---- 

## General data cleaning, manipulation, and visualization
    
    library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

  ## Multivariate data analysis

    library(vegan)

## Loading required package: permute
## Loading required package: lattice
## This is vegan 2.6-4

    library(goeveg)

## This is GoeVeg 0.7.2 - build: 2024-02-06

  ## Multivariate visualization

    library(ecodist)

## 
## Attaching package: 'ecodist'
## 
## The following object is masked from 'package:vegan':
## 
##     mantel

    library(ecotraj)

## Loading required package: Rcpp

  ## Important functions to get the Bray-Curtis dissimilarity code to work

    library(maditr)

## 
## To modify variables or add new variables:
##              let(mtcars, new_var = 42, new_var2 = new_var*hp) %>% head()
## 
## 
## Attaching package: 'maditr'
## 
## The following objects are masked from 'package:dplyr':
## 
##     between, coalesce, first, last
## 
## The following object is masked from 'package:purrr':
## 
##     transpose
## 
## The following object is masked from 'package:readr':
## 
##     cols

  ## Making multipaneled figures

    library(cowplot)

## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp

    library(ggpubr)

## 
## Attaching package: 'ggpubr'
## 
## The following object is masked from 'package:cowplot':
## 
##     get_legend

  ## Needed to draw images in plots

    library(magick)

## Linking to ImageMagick 6.9.12.93
## Enabled features: cairo, fontconfig, freetype, heic, lcms, pango, raw, rsvg, webp
## Disabled features: fftw, ghostscript, x11

  ## Prevents text from overlapping on figures

    library(ggrepel)

  ## Ternary plots

    library(ggtern)

## Registered S3 methods overwritten by 'ggtern':
##   method           from   
##   grid.draw.ggplot ggplot2
##   plot.ggplot      ggplot2
##   print.ggplot     ggplot2
## --
## Remember to cite, run citation(package = 'ggtern') for further info.
## --
## 
## Attaching package: 'ggtern'
## 
## The following objects are masked from 'package:ggplot2':
## 
##     aes, annotate, ggplot, ggplot_build, ggplot_gtable, ggplotGrob,
##     ggsave, layer_data, theme_bw, theme_classic, theme_dark,
##     theme_gray, theme_light, theme_linedraw, theme_minimal, theme_void

### ---- Model fitting and analysis ---- 

# Mixed-models 

  ## Fits mixed models
    library(lme4)

## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack

  ## Functions to analyze significance and R^2 components of mixed-models
    library(lmerTest)

## 
## Attaching package: 'lmerTest'
## 
## The following object is masked from 'package:lme4':
## 
##     lmer
## 
## The following object is masked from 'package:stats':
## 
##     step

  ## Looks at variance components for mixed models

    library(MuMIn)

  ##  Anova function

    library(car)

## Loading required package: carData
## 
## Attaching package: 'car'
## 
## The following object is masked from 'package:dplyr':
## 
##     recode
## 
## The following object is masked from 'package:purrr':
## 
##     some

# Fits negative binomial models

    library(MASS)

## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select

# Post hoc analysis

    library(emmeans)

  ## Below is code to install the pairwiseAdonis package from github

    ### library(devtools)
    ### install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")

    library(pairwiseAdonis)

## Loading required package: cluster

# Indicator species analysis
    
    library(indicspecies)

# Summary stats

  ## 
   library(mosaic)

## Registered S3 method overwritten by 'mosaic':
##   method                           from   
##   fortify.SpatialPolygonsDataFrame ggplot2
## 
## The 'mosaic' package masks several functions from core packages in order to add 
## additional features.  The original behavior of these functions should not be affected by this.
## 
## Attaching package: 'mosaic'
## 
## The following objects are masked from 'package:car':
## 
##     deltaMethod, logit
## 
## The following object is masked from 'package:lmerTest':
## 
##     rand
## 
## The following object is masked from 'package:lme4':
## 
##     factorize
## 
## The following object is masked from 'package:Matrix':
## 
##     mean
## 
## The following object is masked from 'package:cowplot':
## 
##     theme_map
## 
## The following objects are masked from 'package:goeveg':
## 
##     deg2rad, rad2deg
## 
## The following object is masked from 'package:permute':
## 
##     shuffle
## 
## The following objects are masked from 'package:dplyr':
## 
##     count, do, tally
## 
## The following object is masked from 'package:purrr':
## 
##     cross
## 
## The following object is masked from 'package:ggplot2':
## 
##     stat
## 
## The following objects are masked from 'package:stats':
## 
##     binom.test, cor, cor.test, cov, fivenum, IQR, median, prop.test,
##     quantile, sd, t.test, var
## 
## The following objects are masked from 'package:base':
## 
##     max, mean, min, prod, range, sample, sum

Import Data

spring_cover <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Vegetative Cover Data 2019/Cleaned_Spring_Cover_Data_July_2019.xlsx")

fall_cover <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Vegetative Cover Data 2019/Cleaned_Fall_Cover_Data_September_2019.xlsx")


seed_rain <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Seed Rain 2019/Seed_Rain_Data_2019_FINAL.xlsx")

seed_bank <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Seed Bank 2020/Seed_Bank_2020_Cleaned.xlsx")

traits <- read_excel("~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation /Traits and Species Info/SR_SB_VEG_Traits_Dataset.xlsx")

—–

Methods

Data Collection & Experimental Design

Study Sites

Our study sites were at Tucker Prairie Research Area (38○56’53.6” N, 91○59’40.0” W, Callaway County, MO) and at Prairie Fork Conservation Area (38○58’29.7” N, 91○44’03.3” W, Callaway County, MO). Tucker Prairie is a 59-hectare tract of unplowed North American tallgrass claypan prairie. Less than 0.5 % of intact tallgrass prairie ecosystem (i.e., never-been-plowed) remains in Missouri, and Tucker Prairie represents the last sizable remnant prairie in north central Missouri (Samson & Knopf, 1994). More than 250 species of plants inhabit Tucker Prairie, with representatives from 57 families and over 150 genera (Tropicos, n.d.). Native C4 grasses dominate the landscape including Sorghastrum nutans, Andropogon gerardii, Schizachyrium scoparium, and Sporobolus heterolepis (Drew, 1947; Rabinowitz & Rapp, 1980). Before 1957, Tucker Prairie was used for cattle grazing and occasional haying (Drew, 1947). From 1958 to 2002, Tucker Prairie was burned once every four years in the late winter or early spring (Rabinowitz & Rapp, 1980). Since 2002, Tucker Prairie has been managed on a 5-year burn rotation, where units are burned once in the late winter to early spring (Jan. – Mar.) and again 2-3 years later in the late summer to early fall (Aug. – Oct.).

Prairie Fork Conservation Area (PFCA) is over 450 hectares of former agricultural land being restored to tallgrass prairie and savanna ecosystems (Navarrete-Tindall et al., 2007; Newbold et al., 2019). From 2004 onwards, 16-25 hectares are newly seeded each year with native prairie species collected from Tucker Prairie and other nearby remnant prairies (Newbold et al., 2019). As a result, the plots within PFCA represent a chronosequence of reconstructed tallgrass prairies that are mainly comprised of Tucker Prairie descendants. Prior to seeding, glyphosate-resistant crops (i.e., corn and soybeans) are planted and harvested for at least three years to reduce the existing weed seed bank (Newbold et al., 2019). A seed mix containing an average of 179 species of native forbs, legumes, grasses, rushes, and sedges are broadcast seeded at a rate of 13.4 to 18.2 kg/ha once during the dormant season (Newbold et al., 2019). All mixes are comprised of at least 75% of the same species each year that are representative of nearby remnant prairies such as Tucker Prairie (Newbold et al., 2019). Reconstructions are managed using a 2-4-year burn schedule and annual spraying of invasive species (Newbold et al., 2019). For additional details on PFCA management, see Newbold et al. 2019. To capture changes in seed rain dynamics during the restoration process, we grouped our reconstructed sites into three categories, as defined by Newbold et al. (2019) as being compositionally different. We measured the seed rain in an old reconstruction (seeded in 2004), middle aged reconstruction (seeded in 2012 and 2013), and a young reconstruction (seeded in 2017).

Seed rain sampling

In May 2019, we deployed artificial turf grass seed traps (0.1 x 0.1 m) in Tucker Prairie and at each of the focal PFCA restorations. We used artificial turf grass traps instead of sticky traps due to their durability and resistance to freezing (Molau & Per Mølgaard, 1996; Bass, 2015). Turfgrass traps also do not lose effectiveness over time, like sticky traps, which are prone to contamination by soil, non-seed plant debris, and dead insects, making seed recovery difficult (Arruda et al., 2020). Furthermore, unlike funnel traps, they are easy to install and collect with minimal disturbance to surrounding vegetation and soil and do not fill with water (Schott, 1995). At each site, we randomly established ten, 5 m long transects, placed a minimum of 50 m apart from each other. We then placed traps at 1 m intervals along each transect and affixed them to the ground using ground staples (n = 5 traps per transect, 50 traps per site). We collected and replaced seed traps every 2-weeks during the growing season from May 31st to December 12th, 2019. Hereafter we refer to all captured diaspores as “seeds” despite catching a variety of fruit types (e.g., achene). After collection, all traps from a transect and sampling period were grouped together and sieved through a series of soil sieves (1 mm, 500 mm, 250 mm mesh). From these layers, we then picked out and identified seeds to the lowest taxonomic level possible using a reference collection of species inhabiting our study areas and identification guides (Coons et al., n.d.; A. C. Martin & Barkley, 1961).

Because of their extremely small size (0.3 – 0.5 mm; Yatskievych, 1999), we estimated the number of rush seeds (Juncus sp.) when there were over 200 seeds present in a sample. We first sieved the samples through 180 mm and 150 mm mesh soil sieves. Then we subsampled each sieved layer by calculating the average number of rush seeds per 1 cm2 area (n = 3) and multiplying that average by the total area covered by the sample layer. Lastly, we summed the number of estimated rush seeds per layer to calculate the total number of rush seeds in a sample.

Vegetation sampling

In 2019, at the same transects, we sampled species cover twice, once in the early summer (June - July) and again at peak biomass (September). During the first sampling period, we only recorded species presence. At peak biomass, we also sampled percent aerial species cover in a 1m2 area centered around each seed rain trap. We identified species according to Yatskievych (1999).

Seed bank sampling

In March 2020, before the growing season, we collected soil samples to assess the soil seed bank at each focal prairie. To compare the soil seed bank and the previous year’s seed rain, we collected 5 (10 x 10 x 10 cm3; 1000 cm3) soil cores approximately 0.5 m away from where we captured seed rain at each transect in 2019. We allowed the soil samples to dry before removing all non-seed plant material (e.g., roots and rhizomes) to ensure seedlings only germinated from seeds. We then homogenized all samples collected from a particular transect and subsampled ~1500 cm3 of soil to spread ~1 cm deep over sterile potting soil in plastic germination trays. Starting July 2020, we placed the trays (n = 40) in a greenhouse and watered them when dry. Amongst the 40 trays containing soil samples, we also placed an additional 10 trays containing only sterile potting soil to assess whether any contamination occurred from external sources. Once identifiable, we removed seedlings from trays. After germination ceased in July 2021, we placed all trays into a cold room for ~4 months to replicate conditions necessary to break dormancy for any still dormant seeds. Post-vernalization, we returned the trays to the greenhouse, stirred the soil samples, and monitored for any additional germination. We ended the study one year post-vernalization.

—–

General Dataset Cleaning

# Make the six letter species codes for the spring and fall cover uppercase

spring_cover[[3]] <- toupper(spring_cover[[3]])

fall_cover[[4]] <- toupper(fall_cover[[4]])

# Remove unnecessary columns from datasets

### Spring cover - Only keep Site, Transect, SPP6, Scientific Name

spring_cover <- spring_cover[,-c(5,6)]

### Fall cover - Keep everything

### Seed rain - Only keep Plot, Transect, Sampling Date, Week, SPP6, Number, Unknown

seed_rain <- seed_rain[, -c(7,8)]

### Seed bank - Keep everything


# Change names of columns to better reflect datasets (e.g. seed rain: number -> Number_Seeds)

colnames(seed_bank) <- c("Site", "Transect",  "SPP6", "Number_Seedlings")

colnames(seed_rain) <- c("Site", "Transect",  "Sampling_Date", "Week", "SPP6", "Number_Seeds", "Unknown")


# Change transect PFCA 2 - 10A to just 10 for the cover datasets

spring_cover <- spring_cover %>% 
  mutate_if(is.character, 
                 str_replace_all, pattern = "10A", replacement = "10")

fall_cover <- fall_cover %>%
  mutate_if(is.character, 
                 str_replace_all, pattern = "10a", replacement = "10")


# Make Transect and Plot variables factors for all datasets

## Spring cover

spring_cover$Transect <- as.factor(spring_cover$Transect)

## Fall cover

fall_cover$Transect <- as.factor(fall_cover$Transect)
fall_cover$Plot <- as.factor(fall_cover$Plot)

## Seed Bank

seed_bank$Transect <- as.factor(seed_bank$Transect)

## Seed Rain

seed_rain$Transect <- as.factor(seed_rain$Transect)

Lump everything to transect level data (total number of seeds per transect, total number of seedlings per transect, total amount of fall cover per transect)

####### Clean Spring

# 0.) Make a frequency column 

spring_cover$Cover <- rep(1, nrow(spring_cover))

####### Clean Fall

# 0.) Make anything less than 1 = 1

fall_cover <- fall_cover %>% 
  mutate(Cover = ifelse(Cover < 1, 1, Cover))

# 1.) get a sum per species per transect of cover.

fall_sum <- fall_cover %>%
              group_by(Site, Transect, SPP6, Scientific_Name) %>%
              summarize(Cover = sum(Cover))

## `summarise()` has grouped output by 'Site', 'Transect', 'SPP6'. You can
## override using the `.groups` argument.

fall_sum

## # A tibble: 1,215 × 5
## # Groups:   Site, Transect, SPP6 [1,214]
##    Site   Transect SPP6   Scientific_Name          Cover
##    <chr>  <fct>    <chr>  <chr>                    <dbl>
##  1 PFCA 1 1        ACAVIR Acalypha virginica        11.5
##  2 PFCA 1 1        AMBART Ambrosia artemisiifolia    6  
##  3 PFCA 1 1        BARVUL Barbarea vulgaris          3.5
##  4 PFCA 1 1        CHAFAS Chamaecrista fasciculata 388  
##  5 PFCA 1 1        CONCAN Conyza canadensis          1  
##  6 PFCA 1 1        ERISPP Erigeron spp.              2  
##  7 PFCA 1 1        ERIVIL Eriochloa villosa          2  
##  8 PFCA 1 1        EUPSER Eupatorium serotinum       3  
##  9 PFCA 1 1        KUMSTI Kummerowia stipulacea      6  
## 10 PFCA 1 1        LESCAP Lespedeza capitata         1  
## # ℹ 1,205 more rows

####### Seed Rain

# 1.) get a sum per species per transect of cover.

seed_rain_sum <- seed_rain %>%
              group_by(Site, Transect, SPP6) %>%
              summarize(Number_Seeds = sum(Number_Seeds))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

# Clean in the same way you did the seed rain analysis

for(i in 1:nrow(seed_rain_sum)){
  if(seed_rain_sum[i,3] == "AGAFAS"){seed_rain_sum[i,3] <- "AGASPP"}
  if(seed_rain_sum[i,3] == "EUPSER"){seed_rain_sum[i,3] <- "EUPSPP"}
  if(seed_rain_sum[i,3] == "GENPUB"){seed_rain_sum[i,3] <- "GENSPP"}
}

####### Seed Bank

# 1.) get a sum per species per transect of cover.

seed_bank_sum <- seed_bank %>%
              group_by(Site, Transect, SPP6) %>%
              summarize(Number_Seedlings = sum(Number_Seedlings))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

Remove Unknowns from Data sets

####### Remove Spring Unknowns

spring_cover <- filter(spring_cover, !grepl('UNK', SPP6))

# Looks good
list(unique(spring_cover$SPP6))

## [[1]]
##   [1] "ACAVIR" "ACERUB" "ACHMIL" "AGRGIG" "AGRHYE" "ALOCAR" "AMBART" "AMBPSI"
##   [9] "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB" "BAPAUS" "BAPBRA"
##  [17] "BARVUL" "BIDARI" "BLECIL" "BROJAP" "CAMRAD" "CAPBUR" "CARSPP" "CARBIC"
##  [25] "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN" "CORLAN" "CORPAL" "CORSPP"
##  [33] "CORTRI" "CROSAG" "CYPACU" "CYPECH" "DALPUR" "DESILL" "DESSES" "DESSPP"
##  [41] "DICLAN" "DIGISC" "ECHPAL" "ELETEN" "ELYVIR" "ERISPP" "ERYYUC" "EUPCOR"
##  [49] "EUPSER" "EUTGYM" "FESARU" "FESPAR" "FRAVIR" "GALAPA" "GALOBT" "GENPUB"
##  [57] "GERCAR" "HELFLE" "HELHEL" "HELMOL" "HELSPP" "HYPPUN" "JUNBRA" "JUNMAR"
##  [65] "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "LACSER" "LEPVIR" "LESCAP" "LESCUN"
##  [73] "LESVIR" "LIAPYC" "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP"
##  [81] "MELSPP" "MONFIS" "MUHSPP" "MYOVER" "OENBIE" "OXADIL" "PARINT" "PENDIG"
##  [89] "PLAVIR" "POAPRA" "POLSAN" "POLVER" "POTSIM" "PYCPIL" "PYCTEN" "RANABO"
##  [97] "RATPIN" "RHUCOP" "ROSCAR" "ROSSPP" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI"
## [105] "SALAZU" "SCHSCO" "SENMAR" "SETFAB" "SETPAR" "SILANT" "SILINT" "SILLAC"
## [113] "SISSPP" "SOLALT" "SOLCAR" "SOLNEM" "SOLRIG" "SORNUT" "SPHOBT" "SPOHET"
## [121] "STRLEI" "SYMERI" "SYMLAN" "SYMLAT" "SYMNOV" "SYMOBL" "SYMOOL" "SYMORB"
## [129] "SYMPRA" "SYMSPP" "THLARV" "OENSPP" "TRAOHI" "TRIPER" "TRIREP" "TRISPP"
## [137] "ULMSPP" "SCISPP" "VERARV" "VERHEL" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT"
## [145] "ZIZAUR" "SCLTRI"

####### Remove Fall Unknowns

fall_sum <- filter(fall_sum, !grepl('UNK', SPP6))

# Looks good
list(unique(fall_sum$SPP6))

## [[1]]
##   [1] "ACAVIR" "AMBART" "BARVUL" "CHAFAS" "CONCAN" "ERISPP" "ERIVIL" "EUPSER"
##   [9] "KUMSTI" "LESCAP" "PENDIG" "PLAVIR" "RUDHIR" "SETFAB" "SOLALT" "STRLEI"
##  [17] "SYMPIL" "ACHMIL" "CORTRI" "DALPUR" "DIGISC" "EUTGYM" "JUNSPP" "PYCTEN"
##  [25] "RATPIN" "RUDSUB" "SCISPP" "SOLSPE" "SORNUT" "SYMERI" "SYMNOV" "SYMORB"
##  [33] "TRICAM" "BROJAP" "HELMOL" "HYPPUN" "OENBIE" "SOLRIG" "SYMLAN" "AGRHYE"
##  [41] "KUMSTR" "LESCUN" "LIAPYC" "PARINT" "SOLCAR" "TAROFF" "BAPALB" "ERASPE"
##  [49] "JUNVIR" "LESVIR" "SYMPRA" "VERARV" "CERSPP" "ECHCRU" "OXADIL" "PLAMAJ"
##  [57] "SILINT" "SILLAC" "VERSPP" "CORLAN" "CYPECH" "CYPSTR" "EUPCOR" "MEDLUP"
##  [65] "SCHSCO" "BIDARI" "FESPAR" "MONFIS" "RANABO" "SOLNEM" "OENFIL" "ULMPUM"
##  [73] "DESSPP" "ELYCAN" "AGATEN" "BAPBRA" "CARSPP" "CIRALT" "ERYYUC" "GENAND"
##  [81] "KOEMAC" "LACSER" "PYCPIL" "RUMCRI" "SPOHET" "SYMOBL" "VERHEL" "BLECIL"
##  [89] "BOLAST" "CYPACU" "SYMLAT" "SYMOOL" "ECHPAL" "MELSPP" "SYMANO" "SYMLAE"
##  [97] "CORPAL" "ANDGER" "EUPALT" "POAPRA" "ALOCAR" "LIASPP" "LYTALA" "SPOCOM"
## [105] "SPILAC" "PANVIR" "SENMAR" "TRIFLA" "HELSPP" "ELYVIR" "HELHEL" "MYOVER"
## [113] "SALAZU" "TRAOHI" "BAPAUS" "LESPRO" "DESSES" "ROSCAR" "AMOCAN" "LEPDEN"
## [121] "CAMRAD" "RUBSPP" "ACERUB" "ZIZAUR" "DICLAN" "EUPPER" "GALOBT" "LYSLAN"
## [129] "POTSIM" "SETPAR" "SOLMIS" "VIOSAG" "AGAFAS" "ANTNEG" "CROSAG" "LINSUL"
## [137] "RHUGLA" "FRAVIR" "GENPUB" "MUHGLA" "JUNBRA" "ULMSPP" "POLVER" "TRISPP"
## [145] "CORSPP" "POLSAN" "SILANT" "TRIPER" "RHUCOP"

####### Remove Seed Rain Unknowns

# 1.) get a sum per species per transect of cover.


seed_rain_sum <- filter(seed_rain_sum, !grepl('UNK', SPP6))
seed_rain_sum <- filter(seed_rain_sum, !grepl('PANICUM?', SPP6))
seed_rain_sum <- filter(seed_rain_sum, !grepl('NONE', SPP6))

# Looks good
list(unique(seed_rain_sum$SPP6))

## [[1]]
##   [1] "ACAVIR"  "AMBART"  "ANAMIN"  "BOLAST"  "CERSPP"  "CHAFAS"  "CHEALB" 
##   [8] "CONCAN"  "DIGISC"  "ERISPP"  "ERIVIL"  "EUPSPP"  "KUMSPP"  "LEPVIR" 
##  [15] "MEDLUP"  "OXADIL"  "PENDIG"  "PLAVIR"  "PYCTEN"  "RUDHIR"  "SETSPP" 
##  [22] "SOLALT"  "SORNUT"  "SYMSPP"  "THLARV"  "TRIPER"  "VEROSP"  "ALOCAR" 
##  [29] "CORTRI"  "CYPECH"  "DICLAN"  "ERYYUC"  "HYPPUN"  "MONFIS"  "MYOVER" 
##  [36] "OENBIE"  "RATPIN"  "SCIGEO"  "SOLRIG"  "SPHOBT"  "STRLEI"  "TRIFLA" 
##  [43] "ACHMIL"  "BARVUL"  "FESARU"  "GERCAR"  "LESCUN"  "LIAPYC"  "PANCAP" 
##  [50] "SYMNOV"  "BIDARI"  "CARBUS"  "DESSPP"  "ERASPE"  "GALAPA"  "LESVIR" 
##  [57] "POAPRA"  "RUDSUB"  "VERHAS"  "VERSPP"  "ANDGER"  "CARFES"  "ECHCRU" 
##  [64] "EREHIE"  "PLAOCC"  "SILINT"  "AGASPP"  "BROJAP"  "CAPBUR"  "CIRALT" 
##  [71] "CORLAN"  "MELSPP"  "MOLVER"  "SCHSCO"  "SCLTRI"  "CYPSTR"  "FESPAR" 
##  [78] "JUNSPP"  "PYCPIL"  "TAROFF"  "AMATUB"  "LESCAP"  "OENFIL"  "KOEMAC" 
##  [85] "DAUCAR"  "LEUVUL"  "VULOCT"  "CYPACU"  "HORPUS"  "SCIPEN"  "ECHPAL" 
##  [92] "SOLNEM"  "CORPAL"  "LINSUL"  "DIAARM"  "HELMOL"  "PERLON"  "RUMCRI" 
##  [99] "BLECIL"  "LYTALA"  "SPOHET"  "ELESPP1" "SPOCOM"  "CARCEP"  "AGRHYE" 
## [106] "GENSPP"  "HELHEL"  "TRAOHI"  "BAPALB"  "CARANN"  "IPOHED"  "AMOCAN" 
## [113] "SALAZU"  "CARBIC"  "EUPCOR"  "EUTGYM"  "GALOBT"  "RUBSPP"  "SETPAR" 
## [120] "ELESPP"  "LYSLAN"  "POLSPP"  "HYPHIR"  "LYCAME"  "SILANT"  "CARDSP" 
## [127] "CROSAG"  "LOBSPI"  "VIOSAG"

####### Remove Seed Bank Unknowns

seed_bank_sum <- filter(seed_bank_sum, !grepl('UNK', SPP6))

# Looks good
list(sort(unique(seed_bank_sum$SPP6)))

## [[1]]
##  [1] "ABUTHE" "ACAVIR" "ACHMIL" "AGRHYE" "ALOCAR" "AMATUB" "AMBART" "ANDGER"
##  [9] "BARVUL" "CARFES" "CARHIR" "CARPAR" "CERGLO" "CHAFAS" "CIRALT" "CONCAN"
## [17] "CORLAN" "CROSAG" "CYPACU" "CYPSPP" "CYPSTR" "DESSPP" "DICLAN" "DIGISC"
## [25] "DIGSAN" "DIPLAC" "ERASPE" "ERIANN" "ERISTR" "EUPHUM" "EUPSER" "EUTGYM"
## [33] "GERCAR" "HORPUS" "HYPPUN" "IPOHED" "JUNINT" "JUNMAR" "KUMSTI" "KUMSTR"
## [41] "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LYSLAN" "MEDLUP" "MELSPP" "MOLVER"
## [49] "OENSPP" "OXADIL" "PANCAP" "PANDIC" "PENDIG" "PERLON" "PLAOCC" "PLAPUS"
## [57] "PLAVIR" "POACHA" "POAPRA" "POPDEL" "POROLE" "PYCPIL" "PYCTEN" "RATPIN"
## [65] "RUDHIR" "RUMCRI" "SCHSCO" "SETFAB" "SETGLA" "SILANT" "SOLALT" "SOLCAR"
## [73] "SORNUT" "SPHOBT" "STRLEI" "SYMLAE" "SYMPIL" "TAROFF" "THLARV" "TOXRAD"
## [81] "TRAOHI" "TRIFLA" "TRIPER" "TRIREP" "VERHAS" "VERPER" "VERTHA"

—–

Merge Spring and Fall cover

# Full join the two cover data sets
cover_merged <- full_join(spring_cover, fall_sum, by = c("Site", "Transect", "SPP6", "Scientific_Name", "Cover")) %>%
                 mutate_all(~replace(., is.na(.), 0))

## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
##   ..1) ...`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated

# Make sure there is no weirdness in the join

#View(cover_merged) 


### Group by and then select the SPP6 with the max cover between the two datasets

  ### Since Spring only had presence/absence we don't want to inflate the fall cover artificially if we saw that species again

cover_merged_max <- cover_merged   %>%
      group_by(Site, Transect, SPP6, Scientific_Name) %>%
      summarise(Cover=max(Cover))

## `summarise()` has grouped output by 'Site', 'Transect', 'SPP6'. You can
## override using the `.groups` argument.

## Not exactly sure what to do about SYMSPP and KUMSPP (things we couldn't ID in the spring but could in the fall)

    ### Going to drop SYMSPP from the dataset since every plot that had SYMSPP in the spring had identified                    Symphyotrichum in the fall 

cover_merged_max  <- filter(cover_merged_max , !grepl('SYMSPP', SPP6))

    ### Checked KUMSPP only PFCA 2 - 3 and PFCA 2 - 9 did not see Kummerowia again in the fall; will drop KUMSPP for all other plots

          ### Removes all the PFCA 1 KUMSPPs

cover_merged_max <- cover_merged_max[!(cover_merged_max$SPP6%in%"KUMSPP" & cover_merged_max$Site%in% "PFCA 1"),]

         ### Removes all PFCA 2 KUMSPPs that were not in PFCA 2 - 3 and PFCA 2- 9

cover_merged_max <- cover_merged_max[!(cover_merged_max$SPP6%in%"KUMSPP" & cover_merged_max$Site%in% "PFCA 2" & cover_merged_max$Transect%in% c(2,4,5,8)),]


### Check all the species codes

sort(unique(cover_merged_max$SPP6))

##   [1] "ACAVIR" "ACERUB" "ACHMIL" "AGAFAS" "AGATEN" "AGRGIG" "AGRHYE" "ALOCAR"
##   [9] "AMBART" "AMBPSI" "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB"
##  [17] "BAPAUS" "BAPBRA" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP" "CAMRAD"
##  [25] "CAPBUR" "CARBIC" "CARSPP" "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN"
##  [33] "CORLAN" "CORPAL" "CORSPP" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR"
##  [41] "DALPUR" "DESILL" "DESSES" "DESSPP" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL"
##  [49] "ELETEN" "ELYCAN" "ELYVIR" "ERASPE" "ERISPP" "ERIVIL" "ERYYUC" "EUPALT"
##  [57] "EUPCOR" "EUPPER" "EUPSER" "EUTGYM" "FESARU" "FESPAR" "FRAVIR" "GALAPA"
##  [65] "GALOBT" "GENAND" "GENPUB" "GERCAR" "HELFLE" "HELHEL" "HELMOL" "HELSPP"
##  [73] "HYPPUN" "JUNBRA" "JUNMAR" "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "KUMSTI"
##  [81] "KUMSTR" "LACSER" "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LESPRO" "LESVIR"
##  [89] "LIAPYC" "LIASPP" "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP"
##  [97] "MELSPP" "MONFIS" "MUHGLA" "MUHSPP" "MYOVER" "OENBIE" "OENFIL" "OENSPP"
## [105] "OXADIL" "PANVIR" "PARINT" "PENDIG" "PLAMAJ" "PLAVIR" "POAPRA" "POLSAN"
## [113] "POLVER" "POTSIM" "PYCPIL" "PYCTEN" "RANABO" "RATPIN" "RHUCOP" "RHUGLA"
## [121] "ROSCAR" "ROSSPP" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU" "SCHSCO"
## [129] "SCISPP" "SCLTRI" "SENMAR" "SETFAB" "SETPAR" "SILANT" "SILINT" "SILLAC"
## [137] "SISSPP" "SOLALT" "SOLCAR" "SOLMIS" "SOLNEM" "SOLRIG" "SOLSPE" "SORNUT"
## [145] "SPHOBT" "SPILAC" "SPOCOM" "SPOHET" "STRLEI" "SYMANO" "SYMERI" "SYMLAE"
## [153] "SYMLAN" "SYMLAT" "SYMNOV" "SYMOBL" "SYMOOL" "SYMORB" "SYMPIL" "SYMPRA"
## [161] "TAROFF" "THLARV" "TRAOHI" "TRICAM" "TRIFLA" "TRIPER" "TRIREP" "TRISPP"
## [169] "ULMPUM" "ULMSPP" "VERARV" "VERHEL" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT"
## [177] "ZIZAUR"

—–

Merge Cover and Seed Rain

### Cover grouping so that it works with the seed rain dataset

# Code to lump certain species together

cover_to_merge_rain <- cover_merged_max

for(i in 1:nrow(cover_to_merge_rain)){
  
  # Group all Symphyotrichum except SYMNOV

  ## SYMPIL
  ## SYMLAN
  ## SYMLAT
  ## SYMERI
  ## SYMANO
  ## SYMLAE
  ## SYMPRA
  ## SYMOBL
  ## SYMOOL

  
  if(cover_to_merge_rain[i,3] == "SYMPIL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMLAN"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMLAT"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMERI"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMANO"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMLAE"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMPRA"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMOBL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  if(cover_to_merge_rain[i,3] == "SYMOOL"){cover_to_merge_rain[i,3] <- "SYMSPP"}
  

# Group all Carex
  ## CARBIC

   if(cover_to_merge_rain[i,3] == "CARBIC"){cover_to_merge_rain[i,3] <- "CARSPP"} 
  
# Group all Kummerowia
  ## KUMSTI
  ## KUMSTR
  
  if(cover_to_merge_rain[i,3] == "KUMSTI"){cover_to_merge_rain[i,3] <- "KUMSPP"} 
  if(cover_to_merge_rain[i,3] == "KUMSTR"){cover_to_merge_rain[i,3] <- "KUMSPP"} 
  
  
# Group all Juncus
  ## JUNMAR
  ## JUNBBRA
  
  if(cover_to_merge_rain[i,3] == "JUNMAR"){cover_to_merge_rain[i,3] <- "JUNSPP"} 
  if(cover_to_merge_rain[i,3] == "JUNBRA"){cover_to_merge_rain[i,3] <- "JUNSPP"} 
  

# Group all Agalinis
  ## AGATEN
  ## AGAFAS
  
  
  if(cover_to_merge_rain[i,3] == "AGATEN"){cover_to_merge_rain[i,3] <- "AGASPP"} 
  if(cover_to_merge_rain[i,3] == "AGAFAS"){cover_to_merge_rain[i,3] <- "AGASPP"} 
  

# Group all Eleocharis
  ## ELETEN 

  if(cover_to_merge_rain[i,3] == "ELETEN"){cover_to_merge_rain[i,3] <- "ELESPP"} 
  
  
# Group all Gentian
  ## GENPUB
  ## GENAND
  
    if(cover_to_merge_rain[i,3] == "GENPUB"){cover_to_merge_rain[i,3] <- "GENSPP"} 
    if(cover_to_merge_rain[i,3] == "GENAND"){cover_to_merge_rain[i,3] <- "GENSPP"} 
  
# Group all Setaria except Setaria parviflora
  ## SETGLA
  ## SETFAB
  
   if(cover_to_merge_rain[i,3] == "SETGLA"){cover_to_merge_rain[i,3] <- "SETSPP"} 
   if(cover_to_merge_rain[i,3] == "SETFAB"){cover_to_merge_rain[i,3] <- "SETSPP"} 
  
  
# Group all Veronica 
  ## VERARV
  
   if(cover_to_merge_rain[i,3] == "VERARV"){cover_to_merge_rain[i,3] <- "VEROSP"} 
  

# Group all Desmodium
  ## DESSES
  
  if(cover_to_merge_rain[i,3] == "DESSES"){cover_to_merge_rain[i,3] <- "DESSPP"} 
  
  
# Group all Eupatorium

  ## EUPALT
  ## EUPSER
  ## EUPPER
  
if(cover_to_merge_rain[i,3] == "EUPALT"){cover_to_merge_rain[i,3] <- "EUPSPP"} 
if(cover_to_merge_rain[i,3] == "EUPSER"){cover_to_merge_rain[i,3] <- "EUPSPP"} 
if(cover_to_merge_rain[i,3] == "EUPPER"){cover_to_merge_rain[i,3] <- "EUPSPP"} 
  
  
  
# Group all Polygala

  ## POLVER
  ## POLSAN
  
if(cover_to_merge_rain[i,3] == "POLVER"){cover_to_merge_rain[i,3] <- "POLSPP"} 
if(cover_to_merge_rain[i,3] == "POLSAN"){cover_to_merge_rain[i,3] <- "POLSPP"} 

  
}  

### Check that it changed the codes

sort(unique(cover_to_merge_rain$SPP6))

##   [1] "ACAVIR" "ACERUB" "ACHMIL" "AGASPP" "AGRGIG" "AGRHYE" "ALOCAR" "AMBART"
##   [9] "AMBPSI" "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB" "BAPAUS"
##  [17] "BAPBRA" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP" "CAMRAD" "CAPBUR"
##  [25] "CARSPP" "CERSPP" "CHAFAS" "CIRALT" "COMUMB" "CONCAN" "CORLAN" "CORPAL"
##  [33] "CORSPP" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR" "DALPUR" "DESILL"
##  [41] "DESSPP" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL" "ELESPP" "ELYCAN" "ELYVIR"
##  [49] "ERASPE" "ERISPP" "ERIVIL" "ERYYUC" "EUPCOR" "EUPSPP" "EUTGYM" "FESARU"
##  [57] "FESPAR" "FRAVIR" "GALAPA" "GALOBT" "GENSPP" "GERCAR" "HELFLE" "HELHEL"
##  [65] "HELMOL" "HELSPP" "HYPPUN" "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "LACSER"
##  [73] "LEPDEN" "LEPVIR" "LESCAP" "LESCUN" "LESPRO" "LESVIR" "LIAPYC" "LIASPP"
##  [81] "LIASQU" "LINSUL" "LOBSPI" "LYSLAN" "LYTALA" "MEDLUP" "MELSPP" "MONFIS"
##  [89] "MUHGLA" "MUHSPP" "MYOVER" "OENBIE" "OENFIL" "OENSPP" "OXADIL" "PANVIR"
##  [97] "PARINT" "PENDIG" "PLAMAJ" "PLAVIR" "POAPRA" "POLSPP" "POTSIM" "PYCPIL"
## [105] "PYCTEN" "RANABO" "RATPIN" "RHUCOP" "RHUGLA" "ROSCAR" "ROSSPP" "RUBSPP"
## [113] "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU" "SCHSCO" "SCISPP" "SCLTRI" "SENMAR"
## [121] "SETPAR" "SETSPP" "SILANT" "SILINT" "SILLAC" "SISSPP" "SOLALT" "SOLCAR"
## [129] "SOLMIS" "SOLNEM" "SOLRIG" "SOLSPE" "SORNUT" "SPHOBT" "SPILAC" "SPOCOM"
## [137] "SPOHET" "STRLEI" "SYMNOV" "SYMORB" "SYMSPP" "TAROFF" "THLARV" "TRAOHI"
## [145] "TRICAM" "TRIFLA" "TRIPER" "TRIREP" "TRISPP" "ULMPUM" "ULMSPP" "VERHEL"
## [153] "VEROSP" "VERSPP" "VIOSAG" "VIOSPP" "VULOCT" "ZIZAUR"

### Sum everything by species by transect again 


cover_to_merge_rain <- cover_to_merge_rain %>%
              group_by(Site, Transect, SPP6) %>%
              summarize(Cover = sum(Cover))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

### Seed grouping so that it works with the cover dataset

seed_rain_names <- left_join(seed_rain_sum, traits)

## Joining with `by = join_by(SPP6)`

rain_to_merge_cover <- seed_rain_names[,c(1,2,3,4,5)]

for(i in 1:nrow(rain_to_merge_cover)){
  
# Group all Carex
  ## CARBUS
  ## CARBIC
  ## CARANN
  ## CARFES
  ## CARCEP
  
  
  if(rain_to_merge_cover[i,3] == "CARBUS"){rain_to_merge_cover[i,3] <- "CARSPP"}
  if(rain_to_merge_cover[i,3] == "CARBIC"){rain_to_merge_cover[i,3] <- "CARSPP"}
  if(rain_to_merge_cover[i,3] == "CARANN"){rain_to_merge_cover[i,3] <- "CARSPP"}
  if(rain_to_merge_cover[i,3] == "CARFES"){rain_to_merge_cover[i,3] <- "CARSPP"}
  if(rain_to_merge_cover[i,3] == "CARCEP"){rain_to_merge_cover[i,3] <- "CARSPP"}
  
 # Group all Eleocharis
  ## ELESPP1
 
  if(rain_to_merge_cover[i,3] == "ELESPP1"){rain_to_merge_cover[i,3] <- "ELESPP"}
  
  
   # Group all Scirpus

 
  if(rain_to_merge_cover[i,3] == "SCIGEO"){rain_to_merge_cover[i,3] <- "SCISPP"}
  if(rain_to_merge_cover[i,3] == "SCIPEN"){rain_to_merge_cover[i,3] <- "SCISPP"}


}


### Sum everything by species by transect again 


rain_to_merge_cover <- rain_to_merge_cover %>%
              group_by(Site, Transect, SPP6) %>%
              summarize(Number_Seeds = sum(Number_Seeds))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

# Looks good
sort(unique(rain_to_merge_cover$SPP6))

##   [1] "ACAVIR" "ACHMIL" "AGASPP" "AGRHYE" "ALOCAR" "AMATUB" "AMBART" "AMOCAN"
##   [9] "ANAMIN" "ANDGER" "BAPALB" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP"
##  [17] "CAPBUR" "CARDSP" "CARSPP" "CERSPP" "CHAFAS" "CHEALB" "CIRALT" "CONCAN"
##  [25] "CORLAN" "CORPAL" "CORTRI" "CROSAG" "CYPACU" "CYPECH" "CYPSTR" "DAUCAR"
##  [33] "DESSPP" "DIAARM" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL" "ELESPP" "ERASPE"
##  [41] "EREHIE" "ERISPP" "ERIVIL" "ERYYUC" "EUPCOR" "EUPSPP" "EUTGYM" "FESARU"
##  [49] "FESPAR" "GALAPA" "GALOBT" "GENSPP" "GERCAR" "HELHEL" "HELMOL" "HORPUS"
##  [57] "HYPHIR" "HYPPUN" "IPOHED" "JUNSPP" "KOEMAC" "KUMSPP" "LEPVIR" "LESCAP"
##  [65] "LESCUN" "LESVIR" "LEUVUL" "LIAPYC" "LINSUL" "LOBSPI" "LYCAME" "LYSLAN"
##  [73] "LYTALA" "MEDLUP" "MELSPP" "MOLVER" "MONFIS" "MYOVER" "OENBIE" "OENFIL"
##  [81] "OXADIL" "PANCAP" "PENDIG" "PERLON" "PLAOCC" "PLAVIR" "POAPRA" "POLSPP"
##  [89] "PYCPIL" "PYCTEN" "RATPIN" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU"
##  [97] "SCHSCO" "SCISPP" "SCLTRI" "SETPAR" "SETSPP" "SILANT" "SILINT" "SOLALT"
## [105] "SOLNEM" "SOLRIG" "SORNUT" "SPHOBT" "SPOCOM" "SPOHET" "STRLEI" "SYMNOV"
## [113] "SYMSPP" "TAROFF" "THLARV" "TRAOHI" "TRIFLA" "TRIPER" "VERHAS" "VEROSP"
## [121] "VERSPP" "VIOSAG" "VULOCT"

# Drop the scientific names they are causing trouble with the join

cover_to_merge_rain <- cover_to_merge_rain


cover_rain <- full_join(cover_to_merge_rain, rain_to_merge_cover) %>%
                 mutate_all(~replace(., is.na(.), 0))

## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:

# Drop TP transect 7 because of missing data

cover_rain <- cover_rain[!(cover_rain$Site%in%"TP" & cover_rain$Transect%in% "7"),]

—–

Merge Cover and Seed Bank

### Cover grouping so that it works with the seed bank dataset

# Code to lump certain species together

cover_to_merge_bank <- cover_merged_max

for(i in 1:nrow(cover_to_merge_bank)){
  
  # Group all Desmodium together

  ## DESSES

  if(cover_to_merge_bank[i,3] == "DESSES"){cover_to_merge_bank[i,3] <- "DESSPP"}

  # Group all Cyperus together

  ## CYPSTR
  ## CYPECH
  ## CYPACU

  if(cover_to_merge_bank[i,3] == "CYPSTR"){cover_to_merge_bank[i,3] <- "CYPSPP"}
  if(cover_to_merge_bank[i,3] == "CYPECH"){cover_to_merge_bank[i,3] <- "CYPSPP"}
  if(cover_to_merge_bank[i,3] == "CYPACU"){cover_to_merge_bank[i,3] <- "CYPSPP"}

  # Group all Oenothera together

  ## OENFIL
  ## OENBIE
  
  if(cover_to_merge_bank[i,3] == "OENFIL"){cover_to_merge_bank[i,3] <- "OENSPP"}
  if(cover_to_merge_bank[i,3] == "OENBIE"){cover_to_merge_bank[i,3] <- "OENSPP"}



}

cover_to_merge_bank <- cover_to_merge_bank %>% 
  group_by(Site, Transect, SPP6) %>% 
  summarize(Cover = sum(Cover))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

### Seed Bank grouping so that it works with the cover dataset

bank_to_merge_cover <- seed_bank_sum

# Code to lump certain species together

for(i in 1:nrow(bank_to_merge_cover)){
  
  # Group all Cerastium together

  ## CERGLO

  if(bank_to_merge_cover[i,3] == "CERGLO"){bank_to_merge_cover[i,3] <- "CERSPP"}

  # Group all Cyperus together

  ## CYPSTR
  ## CYPACU

  if(bank_to_merge_cover[i,3] == "CYPSTR"){bank_to_merge_cover[i,3] <- "CYPSPP"}
  if(bank_to_merge_cover[i,3] == "CYPACU"){bank_to_merge_cover[i,3] <- "CYPSPP"}
  
  # Group all the Carex together

  if(bank_to_merge_cover[i,3] == "CARFES"){bank_to_merge_cover[i,3] <- "CARSPP"}
  
  # Group all Erigeron together

  ## ERIANN
  ## ERISTR

  if(bank_to_merge_cover[i,3] == "ERIANN"){bank_to_merge_cover[i,3] <- "ERISPP"}
  if(bank_to_merge_cover[i,3] == "ERISTR"){bank_to_merge_cover[i,3] <- "ERISPP"}

  # Change Juncus interior to Juncus sp. (since we were able to ID JUNMAR in the cover but not JUNINT)
  
  if(bank_to_merge_cover[i,3] == "JUNINT"){bank_to_merge_cover[i,3] <- "JUNSPP"}
  
}

bank_to_merge_cover <- bank_to_merge_cover %>% 
  group_by(Site, Transect, SPP6) %>% 
  summarize(Number_Seedlings = sum(Number_Seedlings))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

# Drop the scientific names they are causing trouble with the join



cover_bank <- full_join(cover_to_merge_bank, bank_to_merge_cover) %>%
                 mutate_all(~replace(., is.na(.), 0))

## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:

# Drop TP transect 7 because of missing fall cover data

cover_bank <- cover_bank[!(cover_bank$Site%in%"TP" & cover_bank$Transect%in% "7"),]

—–

Merge Seed Rain and Seed Bank

### Seed Bank grouping so that it works with the cover dataset

bank_to_merge_rain <- seed_bank_sum

# Code to lump certain species together

for(i in 1:nrow(bank_to_merge_rain)){
  
  # Group all Cerastium together

  ## CERGLO

  if(bank_to_merge_rain[i,3] == "CERGLO"){bank_to_merge_rain[i,3] <- "CERSPP"}

  # Group all Cyperus together

  ## CYPSTR
  ## CYPACU
  
  if(bank_to_merge_rain[i,3] == "CYPSTR"){bank_to_merge_rain[i,3] <- "CYPSPP"}
  if(bank_to_merge_rain[i,3] == "CYPACU"){bank_to_merge_rain[i,3] <- "CYPSPP"}
  
  # Group all Erigeron together

  ## ERIANN
  ## ERISTR

  if(bank_to_merge_rain[i,3] == "ERIANN"){bank_to_merge_rain[i,3] <- "ERISPP"}
  if(bank_to_merge_rain[i,3] == "ERISTR"){bank_to_merge_rain[i,3] <- "ERISPP"}

  # Group all the Juncus together
  
  ## JUNINT
  ## JUNMAR
  
  if(bank_to_merge_rain[i,3] == "JUNINT"){bank_to_merge_rain[i,3] <- "JUNSPP"}
  if(bank_to_merge_rain[i,3] == "JUNMAR"){bank_to_merge_rain[i,3] <- "JUNSPP"}
  
  # Group all the Setaria together
  
  if(bank_to_merge_rain[i,3] == "SETFAB"){bank_to_merge_rain[i,3] <- "SETSPP"}
  if(bank_to_merge_rain[i,3] == "SETGLA"){bank_to_merge_rain[i,3] <- "SETSPP"}
  
  # Group all the Symphyotrichum together
  
  ## SYMPIL
  ## SYMLAE
  
  if(bank_to_merge_rain[i,3] == "SYMPIL"){bank_to_merge_rain[i,3] <- "SYMSPP"}
  if(bank_to_merge_rain[i,3] == "SYMLAE"){bank_to_merge_rain[i,3] <- "SYMSPP"}
  if(bank_to_merge_rain[i,3] == "SIMPIL"){bank_to_merge_rain[i,3] <- "SYMSPP"}
  
   # Group all the Kummerowia together
  
  ## KUMSTR
  ## KUMSTI
  
  if(bank_to_merge_rain[i,3] == "KUMSTR"){bank_to_merge_rain[i,3] <- "KUMSPP"}
  if(bank_to_merge_rain[i,3] == "KUMSTI"){bank_to_merge_rain[i,3] <- "KUMSPP"}
  
   # Group all the Veronica together
  
  ## VERPER

  if(bank_to_merge_rain[i,3] == "VERPER"){bank_to_merge_rain[i,3] <- "VEROSP"}

  # Group all the Eupatorium together
  
  ## EUPSER

  if(bank_to_merge_rain[i,3] == "EUPSER"){bank_to_merge_rain[i,3] <- "EUPSPP"}
  
  # Merge all the Cardamine together
  
  ## CARHIR
  ## CARPAR

  if(bank_to_merge_rain[i,3] == "CARHIR"){bank_to_merge_rain[i,3] <- "CARDSP"}
  if(bank_to_merge_rain[i,3] == "CARPAR"){bank_to_merge_rain[i,3] <- "CARDSP"}
}

bank_to_merge_rain <- bank_to_merge_rain %>% 
  group_by(Site, Transect, SPP6) %>% 
  summarize(Number_Seedlings = sum(Number_Seedlings))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

rain_to_merge_bank <- seed_rain_sum

for(i in 1:nrow(rain_to_merge_bank)){

  # Group all Cyperus together

  ## CYPSTR
  ## CYPECH
  ## CYPACU

  if(rain_to_merge_bank[i,3] == "CYPSTR"){rain_to_merge_bank[i,3] <- "CYPSPP"}
  if(rain_to_merge_bank[i,3] == "CYPECH"){rain_to_merge_bank[i,3] <- "CYPSPP"}
  if(rain_to_merge_bank[i,3] == "CYPACU"){rain_to_merge_bank[i,3] <- "CYPSPP"}

  # Group all Oenothera together

  ## OENFIL
  ## OENBIE
  
  if(rain_to_merge_bank[i,3] == "OENFIL"){rain_to_merge_bank[i,3] <- "OENSPP"}
  if(rain_to_merge_bank[i,3] == "OENBIE"){rain_to_merge_bank[i,3] <- "OENSPP"}

}


rain_to_merge_bank <- rain_to_merge_bank %>% 
  group_by(Site, Transect, SPP6) %>% 
  summarize(Number_Seeds = sum(Number_Seeds))

## `summarise()` has grouped output by 'Site', 'Transect'. You can override using
## the `.groups` argument.

### Merge the seed rain and seed bank datasets together
rain_bank <- full_join(bank_to_merge_rain, rain_to_merge_bank) %>%
                 mutate_all(~replace(., is.na(.), 0))

## Joining with `by = join_by(Site, Transect, SPP6)`
## `mutate_all()` ignored the following grouping variables:

—-

Seed Rain and Cover Exploratory Analysis

Wasn’t in the cover at all but we caught seeds

Amaranthus tuberculatus
Cardamine spp.
Chenopodium album
Daucus carota
Dianthus armeria
Erechtites hieracifolia
Hordeum pusillum
Ipomoea hederacea
Leucanthemum vulgare
Mollugo verticillata
Panicum capillare
Persicaria longiseta
Platanus occidentalis
Verbena hastata

Didn’t catch seeds but was in the cover

Notable (> 10% cover observed)

Baptisia australis
Potentilla simplex
Silphium laciniatum
Rosa carolina
Fragraria virginica
Campsis radicans
Parthenium integrifolium
Rhus copallinum
Solanum carolinense
Helianthus spp (the pressed one)
Baptisia bracteata
Antennaria neglecta

Minor ( 1% < cover < 10% observed)

Elymus virginica
Rhus glabra
Trifolium campestre
Verbesina helianthoides
Solidago speciosa
Lespedeza procombens
Liatris spp.
Ranunculus abortivus
Juniperus virginiana
Trifolium spp.
Ulmus spp.
Agrostis gigantea
Lactuca serriola
Symphoricarpos orbiculatus
Ulmus pumila
Acer rubrum
Comandra umbellata
Dalea purpurea
Muhlenbergia spp.
Plantago major
Rosa spp.
Spiranthes lacera
Trifolium repens
Viola spp.
Zizea aurea

Trace (cover 1% observed)

Ambrosia psilostachya
Asclepias syriaca
Desmanthus illinoensis
Elymus canadense
Helenium flexuosum
Lepidium densiflorum
Liatris squarrosa
Muhlenbergia glabrifolia
Oenothera spp. (not biennis)
Panicum virgatum
Senna marilandica
Sisyrinchium spp.
Solidago missouriensis.

Dissimilarity analysis

## # A tibble: 10,959 × 5
## # Groups:   Site, Transect [39]
##    Site   Transect SPP6   Abundance Type       
##    <chr>  <fct>    <chr>      <dbl> <chr>      
##  1 PFCA 1 1        ACAVIR      11.5 Aboveground
##  2 PFCA 1 10       ACAVIR       4   Aboveground
##  3 PFCA 1 2        ACAVIR       4.5 Aboveground
##  4 PFCA 1 3        ACAVIR       2   Aboveground
##  5 PFCA 1 4        ACAVIR       4   Aboveground
##  6 PFCA 1 5        ACAVIR       3   Aboveground
##  7 PFCA 1 6        ACAVIR       2.5 Aboveground
##  8 PFCA 1 7        ACAVIR       3   Aboveground
##  9 PFCA 1 8        ACAVIR       1   Aboveground
## 10 PFCA 1 9        ACAVIR       5   Aboveground
## # ℹ 10,949 more rows

# Code obtained and modified from Lauren and Anu's seed rain vs. vegetation paper


##### run dissimilarity loops
Cover_rain_dis_results <-matrix(nrow=0,ncol=7)
sites <-as.vector(unique(cover_rain_dis$Site))


for(s in 1:length(sites)){
  temp <-subset(cover_rain_dis, Site==sites[s])
  temp$SPP6 <- factor(as.character(temp$SPP6))
  
  transects <-as.vector(unique(temp$Transect))
  
  #for(b in 1:length(blocks)){
    #temp_b<-subset(temp, block==blocks[b])

    #plots<-as.vector(unique(temp_b$plot))
   
      for(t in 1:length(transects)){ 
        
        temp_t <- subset(temp, Transect==transects[t])

        #wide form - with a row for above and belowground
        plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
        
        #relative abundances
        trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"hellinger"))
        names(trans_tot)[1]<-"trt"
        
        trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
        names(trans_pa)[1]<-"trt"
        
        trans_int <- cbind(plot_cast$Type,trunc(plot_cast[,-1])) #integers only
        names(trans_int)[1]<-"trt"
        
        trans_int_norm <- cbind(plot_cast$Type,trunc(decostand(plot_cast[,-1], "normalize")*100)) #integers only and normalize
        names(trans_int)[1]<-"trt"
        
        distance_bray<-vegdist(trans_tot[,-1], method = "bray")
        distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")
        distance_kulczynski<-vegdist(trans_tot[,-1], method = "kulczynski")
        distance_cao<-vegdist(trans_int_norm[,-1], method = "cao")
        distance_morisita<-vegdist(trans_int[,-1], method = "morisita")
        

        new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1],
                   distance_cao[1], distance_morisita[1], distance_kulczynski[1])
        Cover_rain_dis_results <-rbind(Cover_rain_dis_results, new.row)
    }
  print(s/length(sites))
}

## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1

row.names(Cover_rain_dis_results)<-NULL
Cover_rain_dis_results<- data.frame(Cover_rain_dis_results)

names(Cover_rain_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard",
                    "dissimilarity_cao", "dissimilarity_morisita", "dissimilarity_kulczynski")

## 
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = Cover_rain_dis_results)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.148341 -0.033152  0.005915  0.046429  0.104723 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.51363    0.01974  26.018  < 2e-16 ***
## SitePFCA 1  -0.13402    0.02721  -4.925 2.01e-05 ***
## SitePFCA 2  -0.07069    0.02721  -2.598   0.0136 *  
## SitePFCA 3   0.02310    0.02721   0.849   0.4017    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05922 on 35 degrees of freedom
## Multiple R-squared:  0.5505, Adjusted R-squared:  0.512 
## F-statistic: 14.29 on 3 and 35 DF,  p-value: 3.056e-06

## Anova Table (Type III tests)
## 
## Response: dissimilarity_bray
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 2.37434  1  676.96 < 2.2e-16 ***
## Site        0.15036  3   14.29 3.056e-06 ***
## Residuals   0.12276 35                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Cover_rain.mod.bray)
## 
## $Site
##                      diff          lwr          upr     p adj
## PFCA 1-TP     -0.13402090 -0.207406554 -0.060635237 0.0001142
## PFCA 2-TP     -0.07069333 -0.144078991  0.002692326 0.0624794
## PFCA 3-TP      0.02310178 -0.050283878  0.096487439 0.8306569
## PFCA 2-PFCA 1  0.06332756 -0.008100791  0.134755919 0.0975147
## PFCA 3-PFCA 1  0.15712268  0.085694321  0.228551031 0.0000055
## PFCA 3-PFCA 2  0.09379511  0.022366758  0.165223468 0.0060381

## 
## Call:
## lm(formula = dissimilarity_jaccard ~ Site, data = Cover_rain_dis_results)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.104048 -0.037405 -0.006949  0.037294  0.192794 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.54155    0.02038  26.578  < 2e-16 ***
## SitePFCA 1  -0.06193    0.02809  -2.205  0.03412 *  
## SitePFCA 2  -0.07926    0.02809  -2.822  0.00782 ** 
## SitePFCA 3   0.04722    0.02809   1.681  0.10162    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06113 on 35 degrees of freedom
## Multiple R-squared:  0.4353, Adjusted R-squared:  0.3869 
## F-statistic: 8.994 on 3 and 35 DF,  p-value: 0.0001491

## Anova Table (Type III tests)
## 
## Response: dissimilarity_jaccard
##              Sum Sq Df  F value    Pr(>F)    
## (Intercept) 2.63947  1 706.3882 < 2.2e-16 ***
## Site        0.10082  3   8.9943 0.0001491 ***
## Residuals   0.13078 35                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Cover_rain.mod.jac)
## 
## $Site
##                      diff         lwr          upr     p adj
## PFCA 1-TP     -0.06192816 -0.13767379  0.013817463 0.1417750
## PFCA 2-TP     -0.07925589 -0.15500152 -0.003510265 0.0374144
## PFCA 3-TP      0.04721930 -0.02852633  0.122964926 0.3485548
## PFCA 2-PFCA 1 -0.01732773 -0.09105311  0.056397651 0.9203909
## PFCA 3-PFCA 1  0.10914746  0.03542208  0.182872842 0.0017366
## PFCA 3-PFCA 2  0.12647519  0.05274981  0.200200570 0.0002772

Fig. - Bray-Curtis dissimilarity

# Contrasts
standardize <- emmeans::emmeans(Cover_rain.mod.bray, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")

##  contrast    estimate     SE df t.ratio p.value
##  PFCA 1 - TP  -0.1340 0.0272 35  -4.925  0.0001
##  PFCA 2 - TP  -0.0707 0.0272 35  -2.598  0.0371
##  PFCA 3 - TP   0.0231 0.0272 35   0.849  0.7122
## 
## P value adjustment: dunnettx method for 3 tests

## 95% confidence intervals 

emmeans::emmeans(Cover_rain.mod.bray, "Site", type = "response")

##  Site   emmean     SE df lower.CL upper.CL
##  TP      0.514 0.0197 35    0.474    0.554
##  PFCA 1  0.380 0.0187 35    0.342    0.418
##  PFCA 2  0.443 0.0187 35    0.405    0.481
##  PFCA 3  0.537 0.0187 35    0.499    0.575
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.380, 0.443, 0.537, 0.514)
lower <- c(0.342, 0.405,  0.499, 0.474) 
upper <- c(0.418, 0.481, 0.575,  0.554)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Cover_rain_dissimilarity <- data.frame(site, response, lower, upper)

Fig. - Jaccard dissimilarity figure

## 95% confidence intervals 

# Contrasts
standardize <- emmeans::emmeans(Cover_rain.mod.jac, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")

##  contrast    estimate     SE df t.ratio p.value
##  PFCA 1 - TP  -0.0619 0.0281 35  -2.205  0.0892
##  PFCA 2 - TP  -0.0793 0.0281 35  -2.822  0.0216
##  PFCA 3 - TP   0.0472 0.0281 35   1.681  0.2427
## 
## P value adjustment: dunnettx method for 3 tests

emmeans::emmeans(Cover_rain.mod.jac, "Site", type = "response")

##  Site   emmean     SE df lower.CL upper.CL
##  TP      0.542 0.0204 35    0.500    0.583
##  PFCA 1  0.480 0.0193 35    0.440    0.519
##  PFCA 2  0.462 0.0193 35    0.423    0.502
##  PFCA 3  0.589 0.0193 35    0.550    0.628
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.480, 0.462, 0.589, 0.542)
lower <- c(0.440, .423,  0.550, 0.500) 
upper <- c(0.519, 0.502, 0.628,  0.583)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Cover_rain_jac_dissimilarity <- data.frame(site, response, lower, upper)

Immigration

Proportion of new species in the seed rain not seen in local cover

cover_rain_dis_wide <- cover_rain_dis %>% 
  spread(SPP6, Abundance) %>% 
  mutate_all(~replace(., is.na(.), 0))

## `mutate_all()` ignored the following grouping variables:
## • Columns `Site`, `Transect`
## ℹ Use `mutate_at(df, vars(-group_cols()), myoperation)` to silence the message.

cover_rain_dis_long <- cover_rain_dis_wide %>% 
  gather(key = "SPP6", value = "Abundance", 4:177) %>% 
  spread(key = Type, Abundance) 

# Get rid of columns that both have zeros
cover_rain_dis_long_unique <- cover_rain_dis_long[!(cover_rain_dis_long$Aboveground == 0 & cover_rain_dis_long$Seed_Rain == 0),]


# 0 = no seeds/seedlings and 1 = seeds/seedlings
cover_rain_dis_long_unique   <- cover_rain_dis_long_unique   %>% 
 mutate(Aboveground = if_else(Aboveground > 0, 1, 0)) %>% 
 mutate(Seed_Rain = if_else(Seed_Rain > 0, 1, 0))

# Make a new column that adds both the seed rain and aboveground presence/absence columns. Then remove rows that = 2 (shares both species)
cover_rain_dis_long_unique$shared <- cover_rain_dis_long_unique$Aboveground + cover_rain_dis_long_unique$Seed_Rain
  
cover_rain_dis_long_unique2 <- cover_rain_dis_long_unique[!(cover_rain_dis_long_unique$shared == 2),]

cover_rain_dis_long_unique2 <- cover_rain_dis_long_unique2[, -6]


# Count unique taxa per site and transect for the seed bank compared to the seed rain

Unique_species_aboveground_com_rain <- cover_rain_dis_long_unique2[, -5]


Unique_species_aboveground_com_rain <- Unique_species_aboveground_com_rain %>% 
  filter(Aboveground != 0 ) %>% 
  group_by(Site, Transect) %>% 
  summarize(New_Aboveground = length(unique(SPP6)))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

# Count unique taxa per site and transect for the seed rain compared to the seed bank

Unique_species_rain_com_aboveground <- cover_rain_dis_long_unique2[, -4]


Unique_species_rain_com_aboveground <- Unique_species_rain_com_aboveground %>% 
  filter(Seed_Rain != 0 ) %>% 
  group_by(Site, Transect) %>% 
  summarize(New_Rain = length(unique(SPP6)))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

# Count the number of shared taxa between communities

cover_rain_dis_long_shared <- cover_rain_dis_long_unique %>% 
  filter(shared == 2) %>% 
  group_by(Site, Transect) %>% 
  summarize(Shared = length(unique(SPP6)))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

## Merge them together

Species_counts_cover_rain <- full_join(Unique_species_aboveground_com_rain, Unique_species_rain_com_aboveground) %>% 
    mutate_all(~replace(., is.na(.), 0))

## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:

## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
##   ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.

Species_counts_cover_rain  <- full_join(Species_counts_cover_rain, cover_rain_dis_long_shared) %>% 
    mutate_all(~replace(., is.na(.), 0))

## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:

## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
##   ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.

Species_counts_cover_rain$Tot_Richness <- Species_counts_cover_rain$New_Aboveground + Species_counts_cover_rain$New_Rain + Species_counts_cover_rain$Shared


# Count the number of species in the seed rain for each transect

Species_counts_cover_rain$Tot_Rain <- Species_counts_cover_rain$New_Rain + Species_counts_cover_rain$Shared


# Count the number of species in the seed bank for each transect

Species_counts_cover_rain$Tot_Aboveground <- Species_counts_cover_rain$New_Aboveground + Species_counts_cover_rain$Shared


## Immigrants

### What proportion of species dispersing in the seed rain are not found in the local cover?

Species_counts_cover_rain$Prop_Immigrant <- Species_counts_cover_rain$New_Rain / Species_counts_cover_rain$Tot_Rain

favstats(Species_counts_cover_rain$Prop_Immigrant~Species_counts_cover_rain$Site)

##   Species_counts_cover_rain$Site       min        Q1    median        Q3
## 1                         PFCA 1 0.2500000 0.2660970 0.2891738 0.3793860
## 2                         PFCA 2 0.2000000 0.2572674 0.2905844 0.3098639
## 3                         PFCA 3 0.1818182 0.2826705 0.3675214 0.3942669
## 4                             TP 0.2173913 0.2666667 0.3225806 0.3750000
##         max      mean         sd  n missing
## 1 0.5000000 0.3247644 0.08347182 10       0
## 2 0.3333333 0.2826759 0.04123406 10       0
## 3 0.4736842 0.3372457 0.09468692 10       0
## 4 0.4444444 0.3215489 0.07336399  9       0

# No difference in proportion of immigrant species
## *Note* better modeled as a binomial regression?


# Species_counts_cover_rain.binom.mod <- glm(cbind(New_Rain, Tot_Rain) ~ Site, family = binomial, data =Species_counts_cover_rain)
# summary(Species_counts_cover_rain.binom.mod )
# Anova(Species_counts_cover_rain.binom.mod , type = "III")

Species_counts_cover_rain$Site <- factor(Species_counts_cover_rain$Site, levels=c("TP", "PFCA 2", "PFCA 3", "PFCA 1"))

Species_counts_cover_rain.lin.mod <- lm( Prop_Immigrant~ Site, data =Species_counts_cover_rain)
summary(Species_counts_cover_rain.lin.mod)

## 
## Call:
## lm(formula = Prop_Immigrant ~ Site, data = Species_counts_cover_rain)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.155428 -0.054688  0.003038  0.052547  0.175236 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.321549   0.025308  12.705 1.14e-14 ***
## SitePFCA 2  -0.038873   0.034885  -1.114    0.273    
## SitePFCA 3   0.015697   0.034885   0.450    0.656    
## SitePFCA 1   0.003215   0.034885   0.092    0.927    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07592 on 35 degrees of freedom
## Multiple R-squared:  0.07626,    Adjusted R-squared:  -0.002916 
## F-statistic: 0.9632 on 3 and 35 DF,  p-value: 0.421

Anova(Species_counts_cover_rain.lin.mod, type = "III")

## Anova Table (Type III tests)
## 
## Response: Prop_Immigrant
##              Sum Sq Df  F value    Pr(>F)    
## (Intercept) 0.93054  1 161.4255 1.141e-14 ***
## Site        0.01666  3   0.9632     0.421    
## Residuals   0.20176 35                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Contrasts
standardize <- emmeans::emmeans(Species_counts_cover_rain.lin.mod, "Site")
 emmeans::contrast(standardize, "trt.vs.ctrl")

##  contrast    estimate     SE df t.ratio p.value
##  PFCA 2 - TP -0.03887 0.0349 35  -1.114  0.5463
##  PFCA 3 - TP  0.01570 0.0349 35   0.450  0.9162
##  PFCA 1 - TP  0.00322 0.0349 35   0.092  0.9968
## 
## P value adjustment: dunnettx method for 3 tests

## 95% confidence intervals 

emmeans::emmeans(Species_counts_cover_rain.lin.mod, "Site", type = "response")

##  Site   emmean     SE df lower.CL upper.CL
##  TP      0.322 0.0253 35    0.270    0.373
##  PFCA 2  0.283 0.0240 35    0.234    0.331
##  PFCA 3  0.337 0.0240 35    0.289    0.386
##  PFCA 1  0.325 0.0240 35    0.276    0.374
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.325,  0.283, 0.337, 0.322)
lower <- c(0.276, 0.234,  0.289, 0.270) 
upper <- c( 0.374,  0.331, 0.386,  0.373)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Immigration_cover_rain <- data.frame(site, response, lower, upper)

#Re-order the levels
Species_counts_cover_rain$Site <- factor(Species_counts_cover_rain$Site, levels=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"))


prop_new_species_cover_rain.plot <- ggplot()+
  geom_jitter(data=Species_counts_cover_rain, aes(x = Site, y = Prop_Immigrant, color = Site), show.legend = FALSE, width = .2, size = 2.5) +
  geom_point(data = Immigration_cover_rain, aes(x = site, y = response), show.legend = FALSE, size = 3.5)+
  geom_errorbar(data = Immigration_cover_rain, aes(x= site, y = response, ymin = lower, ymax = upper), width = 0.15, cex = 1)+
  theme_classic()+
  labs(x = "", y = "Prop. new species in seed rain")+
  scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
  scale_x_discrete(labels = c ("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant"))+
 theme(text=element_text(size=16), legend.key.size=unit(1, "cm"))+
scale_y_continuous(breaks = seq(0, 1, by = .2), limits = c(0, 1.075))+ 
  geom_vline(xintercept = 3.5, linetype = "longdash")


prop_new_species_cover_rain.plot <- prop_new_species_cover_rain.plot + draw_image(
  "~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/ECHPAL_Sil.png",
  x = .5, y = .5, width = .6)+
  draw_image(
  "~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/SPHOBT_Sil.png",
  x = 1.23, y = .5, width = .6) + 
  annotate("text", x = 1.24, y = .97, label = "Vs.", size = 4)+ 
  annotate("text", x = 4 , y = .97, label = "n.s.", size = 4.5)

prop_new_species_cover_rain.plot

Proportion of seeds in the seed rain not seen in local cover

names(cover_rain_dis_long) <- c("Site", "Transect", "SPP6", "Aboveground_Abundance", "Seed_Rain_Abundance")


cover_rain_dis_long_unique_abundance <- full_join(cover_rain_dis_long_unique, cover_rain_dis_long)

## Joining with `by = join_by(Site, Transect, SPP6)`

cover_rain_dis_long_shared_abundance <- cover_rain_dis_long_unique_abundance %>% 
  filter(shared == 2) %>% 
  group_by(Site, Transect) %>% 
  summarize(Shared_Seeds = sum(Seed_Rain_Abundance))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

cover_rain_dis_long_immigrant_abundance <- cover_rain_dis_long_unique_abundance %>% 
  filter(shared == 1) %>% 
  group_by(Site, Transect) %>% 
  summarize(Immigrant_Seeds = sum(Seed_Rain_Abundance))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

Prop_Immigrant_Seeds <- full_join(cover_rain_dis_long_shared_abundance, cover_rain_dis_long_immigrant_abundance)

## Joining with `by = join_by(Site, Transect)`

Prop_Immigrant_Seeds$Total_Seeds <- Prop_Immigrant_Seeds$Shared_Seeds + Prop_Immigrant_Seeds$Immigrant_Seeds

Prop_Immigrant_Seeds$Prop <-  Prop_Immigrant_Seeds$Immigrant_Seeds / Prop_Immigrant_Seeds$Total_Seeds

Prop_Immigrant_Seeds$Site <- factor(Prop_Immigrant_Seeds$Site, levels=c("TP", "PFCA 1", "PFCA 2", "PFCA 3"))


# Prop_immigrant_cover_rain.binom.mod <- glm(cbind(Immigrant_Seeds, Total_Seeds) ~ Site, family = binomial, data =Prop_Immigrant_Seeds)
# summary(Prop_immigrant_cover_rain.binom.mod)
# Anova(Prop_immigrant_cover_rain.binom.mod, type = "III")

Prop_immigrant_cover_rain.lm.mod <- lm(Prop ~ Site, data = Prop_Immigrant_Seeds)
summary(Prop_immigrant_cover_rain.lm.mod)

## 
## Call:
## lm(formula = Prop ~ Site, data = Prop_Immigrant_Seeds)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.066032 -0.020370 -0.010473  0.005199  0.256420 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.07927    0.01874   4.230  0.00016 ***
## SitePFCA 1  -0.06559    0.02583  -2.539  0.01571 *  
## SitePFCA 2  -0.04865    0.02583  -1.884  0.06795 .  
## SitePFCA 3  -0.04254    0.02583  -1.647  0.10855    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05622 on 35 degrees of freedom
## Multiple R-squared:  0.1644, Adjusted R-squared:  0.09274 
## F-statistic: 2.295 on 3 and 35 DF,  p-value: 0.09484

Anova(Prop_immigrant_cover_rain.lm.mod, type = "3")

## Anova Table (Type III tests)
## 
## Response: Prop
##               Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.056559  1 17.8957 0.0001595 ***
## Site        0.021758  3  2.2948 0.0948394 .  
## Residuals   0.110616 35                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Contrasts
 standardize <- emmeans::emmeans(Prop_immigrant_cover_rain.lm.mod, "Site")
 emmeans::contrast(standardize, "trt.vs.ctrl")

##  contrast    estimate     SE df t.ratio p.value
##  PFCA 1 - TP  -0.0656 0.0258 35  -2.539  0.0426
##  PFCA 2 - TP  -0.0487 0.0258 35  -1.884  0.1691
##  PFCA 3 - TP  -0.0425 0.0258 35  -1.647  0.2572
## 
## P value adjustment: dunnettx method for 3 tests

## 95% confidence intervals 

emmeans::emmeans(Prop_immigrant_cover_rain.lm.mod, "Site", type = "response")

##  Site   emmean     SE df  lower.CL upper.CL
##  TP     0.0793 0.0187 35  0.041231   0.1173
##  PFCA 1 0.0137 0.0178 35 -0.022403   0.0498
##  PFCA 2 0.0306 0.0178 35 -0.005471   0.0667
##  PFCA 3 0.0367 0.0178 35  0.000645   0.0728
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.0137, 0.0306, 0.0367, 0.0793)
lower <- c(-0.022403,  -0.005471, 0.000645, 0.041231) 
upper <- c(0.0498,  0.0667, 0.0728,  0.1173)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Prop_Immigration_cover_rain <- data.frame(site, response, lower, upper)

#Re-order the levels
Prop_Immigrant_Seeds$Site <- factor(Prop_Immigrant_Seeds$Site, levels=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"))


prop_immigrant_seeds_cover_rain.plot <- ggplot()+
  geom_jitter(data=Prop_Immigrant_Seeds, aes(x = Site, y = Prop, color = Site), show.legend = FALSE, width = .2, size = 2.5) +
  theme_classic()+
  geom_point(data = Prop_Immigration_cover_rain, aes(x = site, y = response), show.legend = FALSE, size = 3.5)+
  geom_errorbar(data = Prop_Immigration_cover_rain, aes(x= site, y = response, ymin = lower, ymax = upper), width = 0.15, cex = 1)+
  labs(x = "", y = "Prop. new seeds in seed rain", title = "")+
  scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
  scale_x_discrete(labels = c ("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant"))+
 theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), plot.title = element_text(hjust = 0.5))+
scale_y_continuous(breaks = seq(0, 1, by = .2), limits = c(-0.025, 1.075))+ 
  geom_vline(xintercept = 3.5, linetype = "longdash")


# The two "outlier" points are because we caught a lot of Solidago altissima and Barbarea vulgaris at these points but they weren't in the immediate cover. 

prop_immigrant_seeds_cover_rain.plot <- prop_immigrant_seeds_cover_rain.plot +
  draw_image(
  "~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/ECHPAL_Sil.png",
  x = .5, y = .5, width = .6)+
  draw_image(
  "~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/SPHOBT_Sil.png",
  x = 1.23, y = .5, width = .6)+
  annotate("text", x = 1.24, y = .97, label = "Vs.", size = 4) +
  annotate("text", x = 1 , y = .1, label = "*", size = 9.5)

prop_immigrant_seeds_cover_rain.plot

—-

Seed Bank and Vegetative Cover Exploratory Analysis

 cover_bank_sum <- cover_bank %>% 
   group_by(SPP6) %>% 
   summarize(tot.cover = sum(Cover), tot.seeds = sum(Number_Seedlings))

Germinated but didn’t observed in the cover

Abutilon theophrasti
Amaranthus tuberculatus
Cardamine hirsuta
Cardamine parviflora
Digitaria sanguinalis
Dipsacus laciniatus
Euphorbia humistrata
Hordeum pusillum
Ipomoea hederacea
Mollugo verticillata
Panicum capillare
Panicum dichotomiflorum
Persicaria longiseta
Platanus occidentalis (likely a GH contaminate)
Plantago pusilla
Poa chapmaniana
Populus deltoides (likely a GH contaminate)
Portulaca oleracea
Setaria glauca
Toxicodendron radicans
Verbena hastata
Veronica peregrina
Verbascum thapsus

Almost everything we were able to germinate that was not in the local cover were weedy annual species. Some things like the Cardamine might have been in the cover but have such an early phenology that we missed them. The two tree species might be greenhouse contaminates since the roof was partially open and these trees were right next door.

Observed in the cover but didn’t germinate

Too many to list (112 taxa)

Clearly, only a small fraction of local species survive the transition from cover -> seed rain -> to germinable seed bank.

Few perennial species were present in the germinable seed bank.

Dissimilarity analysis

## # A tibble: 9,984 × 5
## # Groups:   Site, Transect, SPP6 [7,605]
##    Site   Transect SPP6   Type        Abundance
##    <chr>  <fct>    <chr>  <chr>           <dbl>
##  1 PFCA 1 1        ABUTHE Seed_Bank         0  
##  2 PFCA 1 1        ACAVIR Aboveground      11.5
##  3 PFCA 1 1        ACAVIR Seed_Bank         0  
##  4 PFCA 1 1        ACERUB Aboveground       0  
##  5 PFCA 1 1        ACHMIL Aboveground       0  
##  6 PFCA 1 1        ACHMIL Seed_Bank         0  
##  7 PFCA 1 1        AGAFAS Aboveground       0  
##  8 PFCA 1 1        AGATEN Aboveground       0  
##  9 PFCA 1 1        AGRGIG Aboveground       0  
## 10 PFCA 1 1        AGRHYE Aboveground       0  
## # ℹ 9,974 more rows

# Code obtained and modified from Lauren and Anu's seed bank vs. vegetation paper
## Only included Bray-Curtis and Jaccard dissimilarity matrices since the others produced similar results to these two

##### run dissimilarity loops
cover_bank_dis_results <-matrix(nrow=0,ncol=4)
sites <-as.vector(unique(cover_bank_dis$Site))


for(s in 1:length(sites)){
  temp <-subset(cover_bank_dis, Site==sites[s])
  temp$SPP6 <- factor(as.character(temp$SPP6))
  
  transects <-as.vector(unique(temp$Transect))
  
  #for(b in 1:length(blocks)){
    #temp_b<-subset(temp, block==blocks[b])

    #plots<-as.vector(unique(temp_b$plot))
   
      for(t in 1:length(transects)){ 
        
        temp_t <- subset(temp, Transect==transects[t])

        #wide form - with a row for above and belowground
        plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
        
        #relative abundances
        trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"hellinger"))
        names(trans_tot)[1]<-"trt"
        
        trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
        names(trans_pa)[1]<-"trt"
        
        
        distance_bray<-vegdist(trans_tot[,-1], method = "bray")
        distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")

        

        new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1])
        cover_bank_dis_results  <-rbind( cover_bank_dis_results, new.row)
    }
  print(s/length(sites))
}

## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1

row.names(cover_bank_dis_results )<-NULL
cover_bank_dis_results <- data.frame(cover_bank_dis_results )

names(cover_bank_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard")

## 
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = cover_bank_dis_results)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.187831 -0.029607 -0.002413  0.043986  0.171576 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.765696   0.023016  33.268  < 2e-16 ***
## SitePFCA 1  -0.192711   0.031725  -6.074 6.16e-07 ***
## SitePFCA 2  -0.007888   0.031725  -0.249    0.805    
## SitePFCA 3   0.044227   0.031725   1.394    0.172    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06905 on 35 degrees of freedom
## Multiple R-squared:  0.6634, Adjusted R-squared:  0.6345 
## F-statistic: 22.99 on 3 and 35 DF,  p-value: 2.111e-08

## Anova Table (Type III tests)
## 
## Response: dissimilarity_bray
##             Sum Sq Df  F value    Pr(>F)    
## (Intercept) 5.2766  1 1106.754 < 2.2e-16 ***
## Site        0.3289  3   22.993 2.111e-08 ***
## Residuals   0.1669 35                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Cover_bank.mod.bray)
## 
## $Site
##                       diff         lwr         upr     p adj
## PFCA 1-TP     -0.192711214 -0.27827161 -0.10715081 0.0000036
## PFCA 2-TP     -0.007887543 -0.09344794  0.07767286 0.9945027
## PFCA 3-TP      0.044227483 -0.04133292  0.12978788 0.5113778
## PFCA 2-PFCA 1  0.184823671  0.10154529  0.26810205 0.0000047
## PFCA 3-PFCA 1  0.236938697  0.15366032  0.32021707 0.0000000
## PFCA 3-PFCA 2  0.052115027 -0.03116335  0.13539340 0.3452273

## 
## Call:
## lm(formula = dissimilarity_jaccard ~ Site, data = cover_bank_dis_results)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.108134 -0.051282 -0.007532  0.055907  0.133446 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.83530    0.02209  37.818   <2e-16 ***
## SitePFCA 1  -0.05277    0.03045  -1.733   0.0918 .  
## SitePFCA 2  -0.00717    0.03045  -0.236   0.8152    
## SitePFCA 3   0.01687    0.03045   0.554   0.5829    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06626 on 35 degrees of freedom
## Multiple R-squared:  0.1469, Adjusted R-squared:  0.07374 
## F-statistic: 2.008 on 3 and 35 DF,  p-value: 0.1307

## Anova Table (Type III tests)
## 
## Response: dissimilarity_jaccard
##             Sum Sq Df   F value Pr(>F)    
## (Intercept) 6.2796  1 1430.1902 <2e-16 ***
## Site        0.0265  3    2.0083 0.1307    
## Residuals   0.1537 35                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = x)
## 
## $Site
##                       diff         lwr        upr     p adj
## PFCA 1-TP     -0.052771996 -0.13488081 0.02933682 0.3222879
## PFCA 2-TP     -0.007169957 -0.08927877 0.07493886 0.9953158
## PFCA 3-TP      0.016874135 -0.06523468 0.09898295 0.9447732
## PFCA 2-PFCA 1  0.045602039 -0.03431681 0.12552089 0.4259519
## PFCA 3-PFCA 1  0.069646131 -0.01027272 0.14956498 0.1060862
## PFCA 3-PFCA 2  0.024044092 -0.05587476 0.10396294 0.8486613

Fig. - Bray-Curtis dissimilarity

# Contrasts
standardize <- emmeans::emmeans(Cover_bank.mod.bray, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")

##  contrast    estimate     SE df t.ratio p.value
##  PFCA 1 - TP -0.19271 0.0317 35  -6.074  <.0001
##  PFCA 2 - TP -0.00789 0.0317 35  -0.249  0.9752
##  PFCA 3 - TP  0.04423 0.0317 35   1.394  0.3808
## 
## P value adjustment: dunnettx method for 3 tests

## 95% confidence intervals 

emmeans::emmeans(Cover_bank.mod.bray, "Site", type = "response")

##  Site   emmean     SE df lower.CL upper.CL
##  TP      0.766 0.0230 35    0.719    0.812
##  PFCA 1  0.573 0.0218 35    0.529    0.617
##  PFCA 2  0.758 0.0218 35    0.713    0.802
##  PFCA 3  0.810 0.0218 35    0.766    0.854
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.573, 0.758, 0.810, 0.766)
lower <- c(0.529, 0.714,  0.766, 0.719) 
upper <- c(0.617, 0.802, 0.854,  0.812)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Cover_bank_dissimilarity <- data.frame(site, response, lower, upper)

Fig. - Jaccard dissimilarity figure

## 95% confidence intervals 

standardize <- emmeans::emmeans(Cover_bank.mod.jac, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")

##  contrast    estimate     SE df t.ratio p.value
##  PFCA 1 - TP -0.05277 0.0304 35  -1.733  0.2218
##  PFCA 2 - TP -0.00717 0.0304 35  -0.236  0.9778
##  PFCA 3 - TP  0.01687 0.0304 35   0.554  0.8725
## 
## P value adjustment: dunnettx method for 3 tests

emmeans::emmeans(Cover_bank.mod.jac, "Site", type = "response")

##  Site   emmean     SE df lower.CL upper.CL
##  TP      0.835 0.0221 35    0.790    0.880
##  PFCA 1  0.783 0.0210 35    0.740    0.825
##  PFCA 2  0.828 0.0210 35    0.786    0.871
##  PFCA 3  0.852 0.0210 35    0.810    0.895
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.783, 0.828, 0.852, 0.835)
lower <- c(0.740, .786,  0.810, 0.790) 
upper <- c( 0.825, 0.871, 0.895,  0.880)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Cover_bank_jac_dissimilarity <- data.frame(site, response, lower, upper)

Immigrants

Proportion of new species in the seed bank not seen in local cover

cover_bank_dis_wide <- cover_bank_dis %>% 
  spread(SPP6, Abundance) %>% 
  mutate_all(~replace(., is.na(.), 0))

## `mutate_all()` ignored the following grouping variables:
## • Columns `Site`, `Transect`
## ℹ Use `mutate_at(df, vars(-group_cols()), myoperation)` to silence the message.

cover_bank_dis_long <- cover_bank_dis_wide %>% 
  gather(key = "SPP6", value = "Abundance", 4:198) %>% 
  spread(key = Type, Abundance) 

# Get rid of columns that both have zeros
cover_bank_dis_long_unique <- cover_bank_dis_long[!(cover_bank_dis_long$Seed_Bank == 0 & cover_bank_dis_long$Aboveground == 0),]


# 0 = no seeds/seedlings and 1 = seeds/seedlings
cover_bank_dis_long_unique   <- cover_bank_dis_long_unique   %>% 
  mutate(Seed_Bank = if_else(Seed_Bank > 0, 1, 0)) %>% 
   mutate(Aboveground = if_else(Aboveground > 0, 1, 0))

# Make a new column that adds both the cover and seed bank presence/absence columns. Then remove rows that = 2 (shares both species)
cover_bank_dis_long_unique$shared <- cover_bank_dis_long_unique$Seed_Bank + cover_bank_dis_long_unique$Aboveground
  
cover_bank_dis_long_unique2 <- cover_bank_dis_long_unique[!(cover_bank_dis_long_unique$shared == 2),]

cover_bank_dis_long_unique2 <- cover_bank_dis_long_unique2[, -6]


# Count unique taxa per site and transect for the seed bank compared to the cover

Unique_species_bank_com_cover <- cover_bank_dis_long_unique2[, -4]


Unique_species_bank_com_cover <- Unique_species_bank_com_cover %>% 
  filter(Seed_Bank != 0 ) %>% group_by(Site, Transect) %>% 
  summarize(New_Bank = length(unique(SPP6)))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

# Count unique taxa per site and transect for the cover compared to the seed bank

Unique_species_cover_com_bank <- cover_bank_dis_long_unique2[, -5]


Unique_species_cover_com_bank <- Unique_species_cover_com_bank %>% 
  filter(Aboveground != 0 ) %>% 
  group_by(Site, Transect) %>% 
  summarize(New_Cover = length(unique(SPP6)))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

# Count the number of shared taxa between communities

cover_bank_dis_long_shared <- cover_bank_dis_long_unique %>% 
  filter(shared == 2) %>% 
  group_by(Site, Transect) %>% 
  summarize(Shared = length(unique(SPP6)))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

## Merge them together

Species_counts_cover_bank <- full_join(Unique_species_bank_com_cover, Unique_species_cover_com_bank) %>% 
    mutate_all(~replace(., is.na(.), 0))

## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:

## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
##   ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.

Species_counts_cover_bank  <- full_join(Species_counts_cover_bank, cover_bank_dis_long_shared) %>% 
    mutate_all(~replace(., is.na(.), 0))

## Joining with `by = join_by(Site, Transect)`
## `mutate_all()` ignored the following grouping variables:

## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
##   ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.

Species_counts_cover_bank$Tot_Richness <- Species_counts_cover_bank$New_Bank + Species_counts_cover_bank$New_Cover + Species_counts_cover_bank$Shared


# Count the number of species in the seed rain for each transect

Species_counts_cover_bank$Tot_Cover <- Species_counts_cover_bank$New_Cover + Species_counts_cover_bank$Shared


# Count the number of species in the seed bank for each transect

Species_counts_cover_bank$Tot_Bank <- Species_counts_cover_bank$New_Bank + Species_counts_cover_bank$Shared


## Immigrants

Species_counts_cover_bank$Prop_Immigrant <- Species_counts_cover_bank$New_Bank / Species_counts_cover_bank$Tot_Bank

# No difference in proportion of immigrant species

#Re-order the levels
Species_counts_cover_bank$Site <- factor(Species_counts_cover_bank$Site, levels=c("TP", "PFCA 1", "PFCA 2", "PFCA 3"))



Species_counts_cover_bank.lin.mod <- lm( Prop_Immigrant~ Site, data = Species_counts_cover_bank)
summary(Species_counts_cover_bank.lin.mod)

## 
## Call:
## lm(formula = Prop_Immigrant ~ Site, data = Species_counts_cover_bank)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.266163 -0.091067 -0.006404  0.084063  0.283837 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.46616    0.04369  10.669 1.51e-12 ***
## SitePFCA 1  -0.12643    0.06022  -2.099   0.0431 *  
## SitePFCA 2  -0.05650    0.06022  -0.938   0.3546    
## SitePFCA 3   0.10631    0.06022   1.765   0.0863 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1311 on 35 degrees of freedom
## Multiple R-squared:  0.3251, Adjusted R-squared:  0.2673 
## F-statistic: 5.621 on 3 and 35 DF,  p-value: 0.002972

Anova(Species_counts_cover_bank.lin.mod, type = "III")

## Anova Table (Type III tests)
## 
## Response: Prop_Immigrant
##              Sum Sq Df  F value    Pr(>F)    
## (Intercept) 1.95578  1 113.8357 1.515e-12 ***
## Site        0.28971  3   5.6209  0.002972 ** 
## Residuals   0.60132 35                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Contrasts
standardize <- emmeans::emmeans(Species_counts_cover_bank.lin.mod, "Site")
 emmeans::contrast(standardize, "trt.vs.ctrl")

##  contrast    estimate     SE df t.ratio p.value
##  PFCA 1 - TP  -0.1264 0.0602 35  -2.099  0.1110
##  PFCA 2 - TP  -0.0565 0.0602 35  -0.938  0.6572
##  PFCA 3 - TP   0.1063 0.0602 35   1.765  0.2098
## 
## P value adjustment: dunnettx method for 3 tests

## 95% confidence intervals 

emmeans::emmeans(Species_counts_cover_bank.lin.mod, "Site", type = "response")

##  Site   emmean     SE df lower.CL upper.CL
##  TP      0.466 0.0437 35    0.377    0.555
##  PFCA 1  0.340 0.0414 35    0.256    0.424
##  PFCA 2  0.410 0.0414 35    0.326    0.494
##  PFCA 3  0.572 0.0414 35    0.488    0.657
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.340,  0.410, 0.572, 0.466)
lower <- c(0.256, 0.326,  0.488 , 0.377) 
upper <- c( 0.424,  0.494,  0.657,  0.555)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Immigration_cover_bank <- data.frame(site, response, lower, upper)

### the oldest restoration has the greatest proportion of immigrants -> likely caused by a persistent seed bank formed from when this site was farmed prior to restoration. Weedy non-native species likely were displaced from the aboveground vegetation but remained present in the seed bank. This trend is not seen in Tucker Prairie which doesn't have the same past. 

Species_counts_cover_bank$Site <- factor(Species_counts_cover_bank$Site, levels=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"))

prop_new_species_cover_bank.plot <- ggplot() +
  geom_jitter(data = Species_counts_cover_bank, aes(x = Site, y = Prop_Immigrant, color = Site), show.legend = FALSE, width = 0.2, size = 2.5) +
  geom_point(data = Immigration_cover_bank , aes(x = site, y = response), show.legend = FALSE, size = 3.5)+
  geom_errorbar(data = Immigration_cover_bank , aes(x= site, y = response, ymin = lower, ymax = upper), width = 0.15, cex = 1) +
  theme_classic()+
  labs(x = "", y = "Prop. new species in seed bank")+
  scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
  scale_x_discrete(labels = c ("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant"))+
 theme(text=element_text(size=16), legend.key.size=unit(1, "cm"))+
scale_y_continuous(breaks = seq(0, 1, by = .2), limits = c(0, 1.075))+
  geom_vline(xintercept = 3.5, linetype = "longdash")

prop_new_species_cover_bank.plot <- prop_new_species_cover_bank.plot + draw_image(
  "~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/ECHPAL_Sil.png",
  x = .5, y = .5, width = .6)+
  draw_image(
  "~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/AMBART_Sil.png",
  x = 1.15, y = .5, width = .75) + 
  annotate("text", x = 1.24, y = .97, label = "Vs.", size = 4)+ 
  annotate("text", x = 4 , y = .97, label = "n.s.", size = 4.5)

prop_new_species_cover_bank.plot

Proportion of seeds in the seed bank not seen in local cover

names(cover_bank_dis_long) <- c("Site", "Transect", "SPP6", "Aboveground_Abundance", "Seed_Bank_Abundance")


cover_bank_dis_long_unique_abundance <- full_join(cover_bank_dis_long_unique, cover_bank_dis_long)

## Joining with `by = join_by(Site, Transect, SPP6)`

cover_bank_dis_long_shared_abundance <- cover_bank_dis_long_unique_abundance %>% 
  filter(shared == 2) %>% 
  group_by(Site, Transect) %>% 
  summarize(Shared_Seedlings = sum(Seed_Bank_Abundance))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

cover_bank_dis_long_immigrant_abundance <- cover_bank_dis_long_unique_abundance %>% 
  filter(shared == 1) %>% 
  group_by(Site, Transect) %>% 
  summarize(Immigrant_Seedlings = sum(Seed_Bank_Abundance))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

Prop_Immigrant_Seedlings <- full_join(cover_bank_dis_long_shared_abundance, cover_bank_dis_long_immigrant_abundance)

## Joining with `by = join_by(Site, Transect)`

Prop_Immigrant_Seedlings$Total_Seedlings <- Prop_Immigrant_Seedlings$Shared_Seedlings + Prop_Immigrant_Seedlings$Immigrant_Seedlings

Prop_Immigrant_Seedlings$Prop <-  Prop_Immigrant_Seedlings$Immigrant_Seedlings / Prop_Immigrant_Seedlings$Total_Seedlings

# No difference in proportion of immigrant species
## *Note* better modeled as a binomial regression?
## binomial regression produces results that don't work with the actual observed data. Linear model works just as well with more appropriate mean and CI results.


#Re-order the levels
Prop_Immigrant_Seedlings$Site <- factor(Prop_Immigrant_Seedlings$Site, levels=c("TP", "PFCA 1", "PFCA 2", "PFCA 3"))



Prop_Immigrant_Seedlings.lin.mod <- lm(Prop ~ Site, data = Prop_Immigrant_Seedlings)
summary(Prop_Immigrant_Seedlings.lin.mod)

## 
## Call:
## lm(formula = Prop ~ Site, data = Prop_Immigrant_Seedlings)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.33648 -0.14610 -0.00805  0.08538  0.34960 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.45040    0.05765   7.812 3.53e-09 ***
## SitePFCA 1  -0.39907    0.07947  -5.022 1.50e-05 ***
## SitePFCA 2  -0.19101    0.07947  -2.404   0.0217 *  
## SitePFCA 3   0.11288    0.07947   1.420   0.1643    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.173 on 35 degrees of freedom
## Multiple R-squared:  0.5891, Adjusted R-squared:  0.5539 
## F-statistic: 16.72 on 3 and 35 DF,  p-value: 6.561e-07

Anova(Prop_Immigrant_Seedlings.lin.mod, type = "III")

## Anova Table (Type III tests)
## 
## Response: Prop
##             Sum Sq Df F value    Pr(>F)    
## (Intercept) 1.8258  1  61.032 3.529e-09 ***
## Site        1.5009  3  16.725 6.561e-07 ***
## Residuals   1.0470 35                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Contrasts
standardize <- emmeans::emmeans(Prop_Immigrant_Seedlings.lin.mod , "Site")
 emmeans::contrast(standardize, "trt.vs.ctrl")

##  contrast    estimate     SE df t.ratio p.value
##  PFCA 1 - TP   -0.399 0.0795 35  -5.022  <.0001
##  PFCA 2 - TP   -0.191 0.0795 35  -2.404  0.0579
##  PFCA 3 - TP    0.113 0.0795 35   1.420  0.3666
## 
## P value adjustment: dunnettx method for 3 tests

## 95% confidence intervals 

emmeans::emmeans(Prop_Immigrant_Seedlings.lin.mod, "Site", type = "response")

##  Site   emmean     SE df lower.CL upper.CL
##  TP     0.4504 0.0577 35   0.3334    0.567
##  PFCA 1 0.0513 0.0547 35  -0.0597    0.162
##  PFCA 2 0.2594 0.0547 35   0.1484    0.370
##  PFCA 3 0.5633 0.0547 35   0.4522    0.674
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.0513,  0.2594, 0.5633, 0.4504)
lower <- c(-0.0597, 0.1484,  0.4522 , 0.3334) 
upper <- c( 0.162,  0.370,  0.674,  0.567)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Prop_cover_bank <- data.frame(site, response, lower, upper)

#Re-order the levels
Prop_Immigrant_Seedlings$Site <- factor(Prop_Immigrant_Seedlings$Site, levels=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"))


prop_immigrant_seedlings_cover_bank.plot <- ggplot()+
  geom_jitter(data=Prop_Immigrant_Seedlings, aes(x = Site, y = Prop, color = Site), show.legend = FALSE, width = .2, size = 2.5) +
  geom_point(data = Prop_cover_bank, aes(x = site, y = response), show.legend = FALSE, size = 3.5)+
  geom_errorbar(data = Prop_cover_bank, aes(x= site, y = response, ymin = lower, ymax = upper), width = 0.15, cex = 1) +
  theme_classic()+
  labs(x = "", y = "Prop. new seeds in seed bank", title = "")+
  scale_color_manual(name = "Site", values=c("#E69F00", "#009E73", "#56B4E9", "#0072B2"), breaks=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"), labels = c("Young", "Middle", "Old", "Remnant"))+
  scale_x_discrete(labels = c ("PFCA 1" = "Young", "PFCA 2" = "Middle", "PFCA 3" = "Old", "TP" = "Remnant"))+
 theme(text=element_text(size=16), legend.key.size=unit(1, "cm"), plot.title = element_text(hjust = 0.5))+
scale_y_continuous(breaks = seq(0, 1, by = .2), limits = c(-0.06, 1.075))+
  geom_vline(xintercept = 3.5, linetype = "longdash")



prop_immigrant_seedlings_cover_bank.plot <- prop_immigrant_seedlings_cover_bank.plot + draw_image(
  "~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/ECHPAL_Sil.png",
  x = .5, y = .5, width = .6)+
  draw_image(
  "~/SR_SB_Veg_Analysis/Chapter_2.Seed_Bank_Seed Rain_and_Aboveground_Vegetation/Images/AMBART_Sil.png",
  x = 1.15, y = .5, width = .75) + 
  annotate("text", x = 1.24, y = .97, label = "Vs.", size = 4)+ 
  annotate("text", x = 1 , y = .4, label = "*", size = 9.5)

—-

Seed Rain and Seed Bank Exploratory Analysis

rain_bank_sum <- rain_bank %>% 
  group_by(SPP6) %>% 
  summarize(tot.seeds = sum(Number_Seeds), tot.seedlings = sum(Number_Seedlings))

Germinated but didn’t catch in the seed rain

Abutilon theophrasti
Digitalis sanguinalis
Dipsacus laciniatus
Euphorbia humistrata
Lepidium densiflorum
Panicum dichotomiflorum
Plantago pusilla
Poa chapmaniana

Mostly ruderal annual species.

Caught in the seed rain but didn’t germinate

Too many to list (61 taxa)

Dissimilarity analysis

## # A tibble: 9,984 × 5
## # Groups:   Site, Transect, SPP6 [7,605]
##    Site   Transect SPP6   Type        Abundance
##    <chr>  <fct>    <chr>  <chr>           <dbl>
##  1 PFCA 1 1        ABUTHE Seed_Bank         0  
##  2 PFCA 1 1        ACAVIR Aboveground      11.5
##  3 PFCA 1 1        ACAVIR Seed_Bank         0  
##  4 PFCA 1 1        ACERUB Aboveground       0  
##  5 PFCA 1 1        ACHMIL Aboveground       0  
##  6 PFCA 1 1        ACHMIL Seed_Bank         0  
##  7 PFCA 1 1        AGAFAS Aboveground       0  
##  8 PFCA 1 1        AGATEN Aboveground       0  
##  9 PFCA 1 1        AGRGIG Aboveground       0  
## 10 PFCA 1 1        AGRHYE Aboveground       0  
## # ℹ 9,974 more rows

## # A tibble: 8,200 × 5
## # Groups:   Site, Transect, SPP6 [5,600]
##    Site   Transect SPP6   Type      Abundance
##    <chr>  <fct>    <chr>  <chr>         <dbl>
##  1 PFCA 1 1        ABUTHE Seed_Bank         0
##  2 PFCA 1 1        ACAVIR Seed_Bank         0
##  3 PFCA 1 1        ACAVIR Seed_Rain        22
##  4 PFCA 1 1        ACHMIL Seed_Bank         0
##  5 PFCA 1 1        ACHMIL Seed_Rain         0
##  6 PFCA 1 1        AGASPP Seed_Rain         0
##  7 PFCA 1 1        AGRHYE Seed_Bank         0
##  8 PFCA 1 1        AGRHYE Seed_Rain         0
##  9 PFCA 1 1        ALOCAR Seed_Bank         0
## 10 PFCA 1 1        ALOCAR Seed_Rain         0
## # ℹ 8,190 more rows

# Code obtained and modified from Lauren and Anu's seed bank vs. vegetation paper
## Only included Bray-Curtis and Jaccard dissimilarity matrices since the others produced similar results to these two

##### run dissimilarity loops
rain_bank_dis_results <-matrix(nrow=0,ncol=4)
sites <-as.vector(unique(rain_bank_dis$Site))


for(s in 1:length(sites)){
  temp <-subset(rain_bank_dis, Site==sites[s])
  temp$SPP6 <- factor(as.character(temp$SPP6))
  
  transects <-as.vector(unique(temp$Transect))
  
for(t in 1:length(transects)){ 
      
        temp_t <- subset(temp, Transect==transects[t])

        #wide form - with a row for above and belowground
        plot_cast<-dcast(temp_t[,-c(1:2)],Type ~ SPP6, value.var="Abundance", fun.aggregate=mean,fill=0)
        
        #relative abundances
        trans_tot <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"hellinger"))
        names(trans_tot)[1]<-"trt"
        
        trans_pa <- cbind(plot_cast$Type,decostand(plot_cast[,-1],"pa"))
        names(trans_pa)[1]<-"trt"
        
        
        distance_bray<-vegdist(trans_tot[,-1], method = "bray")
        distance_jaccard <- vegdist(trans_pa[,-1], method = "jaccard")

        

        new.row<-c(sites[s], transects[t],distance_bray[1],distance_jaccard[1])
        rain_bank_dis_results  <-rbind(rain_bank_dis_results, new.row)
    }
  print(s/length(sites))
}

## [1] 0.25
## [1] 0.5
## [1] 0.75
## [1] 1

row.names(rain_bank_dis_results)<-NULL
rain_bank_dis_results <- data.frame(rain_bank_dis_results)

names(rain_bank_dis_results) <- c("Site", "Transect", "dissimilarity_bray", "dissimilarity_jaccard")

## 
## Call:
## lm(formula = dissimilarity_bray ~ Site, data = rain_bank_dis_results)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.170306 -0.047042 -0.002914  0.046581  0.125648 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.61439    0.02451  25.063   <2e-16 ***
## SitePFCA 1  -0.11376    0.03467  -3.281   0.0023 ** 
## SitePFCA 2   0.05321    0.03467   1.535   0.1336    
## SitePFCA 3   0.06543    0.03467   1.887   0.0672 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07752 on 36 degrees of freedom
## Multiple R-squared:  0.481,  Adjusted R-squared:  0.4377 
## F-statistic: 11.12 on 3 and 36 DF,  p-value: 2.602e-05

## Anova Table (Type III tests)
## 
## Response: dissimilarity_bray
##             Sum Sq Df F value    Pr(>F)    
## (Intercept) 3.7748  1  628.18 < 2.2e-16 ***
## Site        0.2005  3   11.12 2.602e-05 ***
## Residuals   0.2163 36                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = Rain_bank.mod.bray)
## 
## $Site
##                      diff         lwr         upr     p adj
## PFCA 1-TP     -0.11375660 -0.20712329 -0.02038992 0.0117875
## PFCA 2-TP      0.05320739 -0.04015930  0.14657407 0.4279537
## PFCA 3-TP      0.06543023 -0.02793646  0.15879691 0.2513341
## PFCA 2-PFCA 1  0.16696399  0.07359731  0.26033068 0.0001492
## PFCA 3-PFCA 1  0.17918683  0.08582015  0.27255352 0.0000514
## PFCA 3-PFCA 2  0.01222284 -0.08114385  0.10558953 0.9847184

Fig. - Bray-Curtis dissimilarity

## 95% confidence intervals 

# Contrasts
standardize <- emmeans::emmeans(Rain_bank.mod.bray, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")

##  contrast    estimate     SE df t.ratio p.value
##  PFCA 1 - TP  -0.1138 0.0347 36  -3.281  0.0065
##  PFCA 2 - TP   0.0532 0.0347 36   1.535  0.3079
##  PFCA 3 - TP   0.0654 0.0347 36   1.887  0.1674
## 
## P value adjustment: dunnettx method for 3 tests

emmeans::emmeans(Rain_bank.mod.bray, "Site", type = "response")

##  Site   emmean     SE df lower.CL upper.CL
##  TP      0.614 0.0245 36    0.565    0.664
##  PFCA 1  0.501 0.0245 36    0.451    0.550
##  PFCA 2  0.668 0.0245 36    0.618    0.717
##  PFCA 3  0.680 0.0245 36    0.630    0.730
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.501, 0.668, 0.680, 0.614)
lower <- c(0.451, 0.618,  0.630, 0.565) 
upper <- c(0.550, 0.717, 0.730,  0.664)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Rain_bank_dissimilarity <- data.frame(site, response, lower, upper)

Fig. - Jaccard dissimilarity figure

#Re-order the levels
rain_bank_dis_results$Site <- factor(rain_bank_dis_results$Site, levels=c("TP", "PFCA 1", "PFCA 2", "PFCA 3"))



Rain_bank.mod.jac <- lm(dissimilarity_jaccard ~ Site, data= rain_bank_dis_results)
summary(Rain_bank.mod.jac)

## 
## Call:
## lm(formula = dissimilarity_jaccard ~ Site, data = rain_bank_dis_results)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.16845 -0.03871  0.00568  0.02850  0.13258 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.76505    0.02153  35.540   <2e-16 ***
## SitePFCA 1  -0.04397    0.03044  -1.444    0.157    
## SitePFCA 2  -0.01655    0.03044  -0.544    0.590    
## SitePFCA 3   0.04111    0.03044   1.350    0.185    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.06807 on 36 degrees of freedom
## Multiple R-squared:  0.1856, Adjusted R-squared:  0.1178 
## F-statistic: 2.735 on 3 and 36 DF,  p-value: 0.05776

Anova(Rain_bank.mod.jac, type = "III")

## Anova Table (Type III tests)
## 
## Response: dissimilarity_jaccard
##             Sum Sq Df   F value  Pr(>F)    
## (Intercept) 5.8529  1 1263.0564 < 2e-16 ***
## Site        0.0380  3    2.7354 0.05776 .  
## Residuals   0.1668 36                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## 95% confidence intervals 

# Contrasts
standardize <- emmeans::emmeans(Rain_bank.mod.jac, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")

##  contrast    estimate     SE df t.ratio p.value
##  PFCA 1 - TP  -0.0440 0.0304 36  -1.444  0.3536
##  PFCA 2 - TP  -0.0166 0.0304 36  -0.544  0.8772
##  PFCA 3 - TP   0.0411 0.0304 36   1.350  0.4047
## 
## P value adjustment: dunnettx method for 3 tests

emmeans::emmeans(Rain_bank.mod.jac, "Site", type = "response")

##  Site   emmean     SE df lower.CL upper.CL
##  TP      0.765 0.0215 36    0.721    0.809
##  PFCA 1  0.721 0.0215 36    0.677    0.765
##  PFCA 2  0.748 0.0215 36    0.705    0.792
##  PFCA 3  0.806 0.0215 36    0.762    0.850
## 
## Confidence level used: 0.95

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c(0.721, 0.748, 0.806, 0.765)
lower <- c(0.677, 0.705,  0.762, 0.721) 
upper <- c(0.765, 0.792, 0.850,  0.809)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")


Rain_bank_jac_dissimilarity <- data.frame(site, response, lower, upper)

Immigrants

Proportion of new species in the seed bank not seen in the seed rain

## Warning: There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
##   ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.
## There were 4 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Transect = (structure(function (..., .x = ..1, .y = ..2, . =
##   ..1) ...`.
## ℹ In group 1: `Site = "PFCA 1"`.
## Caused by warning in `[<-.factor`:
## ! invalid factor level, NA generated
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 3 remaining warnings.

## 
## Call:
## lm(formula = Prop_Immigrant ~ Site, data = Species_counts_rain_bank)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.288874 -0.079552 -0.009654  0.067853  0.220249 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.28887    0.03892   7.423 9.23e-09 ***
## SitePFCA 1  -0.07665    0.05504  -1.393 0.172269    
## SitePFCA 2  -0.06449    0.05504  -1.172 0.248973    
## SitePFCA 3   0.22421    0.05504   4.074 0.000243 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1231 on 36 degrees of freedom
## Multiple R-squared:  0.5179, Adjusted R-squared:  0.4778 
## F-statistic: 12.89 on 3 and 36 DF,  p-value: 7.119e-06

## Anova Table (Type III tests)
## 
## Response: Prop_Immigrant
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.83448  1  55.098 9.227e-09 ***
## Site        0.58579  3  12.893 7.119e-06 ***
## Residuals   0.54524 36                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##  contrast    estimate    SE df t.ratio p.value
##  PFCA 1 - TP  -0.0766 0.055 36  -1.393  0.3812
##  PFCA 2 - TP  -0.0645 0.055 36  -1.172  0.5104
##  PFCA 3 - TP   0.2242 0.055 36   4.074  0.0007
## 
## P value adjustment: dunnettx method for 3 tests

##  Site   emmean     SE df lower.CL upper.CL
##  TP      0.289 0.0389 36    0.210    0.368
##  PFCA 1  0.212 0.0389 36    0.133    0.291
##  PFCA 2  0.224 0.0389 36    0.145    0.303
##  PFCA 3  0.513 0.0389 36    0.434    0.592
## 
## Confidence level used: 0.95

Proportion of seedlings in the seed bank not seen in seed rain

## 
## Call:
## lm(formula = Prop ~ Site, data = Prop_Immigrant_Seedlings_Rain)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32186 -0.07339 -0.02271  0.03860  0.47470 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.28530    0.05768   4.946 1.77e-05 ***
## SitePFCA 1  -0.25212    0.08157  -3.091  0.00384 ** 
## SitePFCA 2  -0.18605    0.08157  -2.281  0.02857 *  
## SitePFCA 3   0.17718    0.08157   2.172  0.03651 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1824 on 36 degrees of freedom
## Multiple R-squared:  0.4845, Adjusted R-squared:  0.4415 
## F-statistic: 11.28 on 3 and 36 DF,  p-value: 2.312e-05

## Anova Table (Type III tests)
## 
## Response: Prop
##              Sum Sq Df F value    Pr(>F)    
## (Intercept) 0.81395  1  24.466 1.770e-05 ***
## Site        1.12545  3  11.277 2.312e-05 ***
## Residuals   1.19766 36                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##  contrast    estimate     SE df t.ratio p.value
##  PFCA 1 - TP   -0.252 0.0816 36  -3.091  0.0108
##  PFCA 2 - TP   -0.186 0.0816 36  -2.281  0.0754
##  PFCA 3 - TP    0.177 0.0816 36   2.172  0.0951
## 
## P value adjustment: dunnettx method for 3 tests

##  Site   emmean     SE df lower.CL upper.CL
##  TP     0.2853 0.0577 36   0.1683    0.402
##  PFCA 1 0.0332 0.0577 36  -0.0838    0.150
##  PFCA 2 0.0992 0.0577 36  -0.0177    0.216
##  PFCA 3 0.4625 0.0577 36   0.3455    0.579
## 
## Confidence level used: 0.95

—–

All Communities

##   [1] "ACAVIR" "ACERUB" "ACHMIL" "AGASPP" "AGRGIG" "AGRHYE" "ALOCAR" "AMATUB"
##   [9] "AMBART" "AMBPSI" "AMOCAN" "ANAMIN" "ANDGER" "ANTNEG" "ASCSYR" "BAPALB"
##  [17] "BAPAUS" "BAPBRA" "BARVUL" "BIDARI" "BLECIL" "BOLAST" "BROJAP" "CAMRAD"
##  [25] "CAPBUR" "CARDSP" "CARSPP" "CERSPP" "CHAFAS" "CHEALB" "CIRALT" "COMUMB"
##  [33] "CONCAN" "CORLAN" "CORPAL" "CORSPP" "CORTRI" "CROSAG" "CYPSPP" "DALPUR"
##  [41] "DAUCAR" "DESILL" "DESSPP" "DIAARM" "DICLAN" "DIGISC" "ECHCRU" "ECHPAL"
##  [49] "ELESPP" "ELYCAN" "ELYVIR" "ERASPE" "EREHIE" "ERISPP" "ERIVIL" "ERYYUC"
##  [57] "EUPCOR" "EUPSPP" "EUTGYM" "FESARU" "FESPAR" "FRAVIR" "GALAPA" "GALOBT"
##  [65] "GENSPP" "GERCAR" "HELFLE" "HELHEL" "HELMOL" "HELSPP" "HORPUS" "HYPHIR"
##  [73] "HYPPUN" "IPOHED" "JUNSPP" "JUNVIR" "KOEMAC" "KUMSPP" "LACSER" "LEPDEN"
##  [81] "LEPVIR" "LESCAP" "LESCUN" "LESPRO" "LESVIR" "LEUVUL" "LIAPYC" "LIASPP"
##  [89] "LIASQU" "LINSUL" "LOBSPI" "LYCAME" "LYSLAN" "LYTALA" "MEDLUP" "MELSPP"
##  [97] "MOLVER" "MONFIS" "MUHGLA" "MUHSPP" "MYOVER" "OENSPP" "OXADIL" "PANCAP"
## [105] "PANVIR" "PARINT" "PENDIG" "PERLON" "PLAMAJ" "PLAOCC" "PLAVIR" "POAPRA"
## [113] "POLSPP" "POTSIM" "PYCPIL" "PYCTEN" "RANABO" "RATPIN" "RHUCOP" "RHUGLA"
## [121] "ROSCAR" "ROSSPP" "RUBSPP" "RUDHIR" "RUDSUB" "RUMCRI" "SALAZU" "SCHSCO"
## [129] "SCISPP" "SCLTRI" "SENMAR" "SETPAR" "SETSPP" "SILANT" "SILINT" "SILLAC"
## [137] "SISSPP" "SOLALT" "SOLCAR" "SOLMIS" "SOLNEM" "SOLRIG" "SOLSPE" "SORNUT"
## [145] "SPHOBT" "SPILAC" "SPOCOM" "SPOHET" "STRLEI" "SYMNOV" "SYMORB" "SYMSPP"
## [153] "TAROFF" "THLARV" "TRAOHI" "TRICAM" "TRIFLA" "TRIPER" "TRIREP" "TRISPP"
## [161] "ULMPUM" "ULMSPP" "VERHAS" "VERHEL" "VEROSP" "VERSPP" "VIOSAG" "VIOSPP"
## [169] "VULOCT" "ZIZAUR" "PLAPUS" "POACHA" "ABUTHE" "POPDEL" "EUPHUM" "PANDIC"
## [177] "POROLE" "VERTHA" "DIPLAC" "DIGSAN" "TOXRAD"

Seedling quantity

### Number of seedlings at each site
Site_seedlings <- everything_merged %>% 
  filter(Type == "Seed_Bank") %>% 
  group_by(Site) %>% 
  summarize(tot.seedlings = sum(Abundance))

### Number of seedlings at each transect and site
Just_seedlings <- everything_merged %>% 
  filter(Type == "Seed_Bank") %>% 
  group_by(Site, Transect) %>% 
  summarize(tot.seedlings = sum(Abundance))

## `summarise()` has grouped output by 'Site'. You can override using the
## `.groups` argument.

## Mean number and standard deviation of seedlings
Just_seedlings_mean <- Just_seedlings %>% 
  group_by(Site) %>% 
  summarize(mean.seedlings = mean(tot.seedlings), sd.seedlings = sd(tot.seedlings))


#Re-order the levels
Just_seedlings$Site <- factor(Just_seedlings$Site, levels=c("TP", "PFCA 1", "PFCA 2", "PFCA 3"))



# Fit a negative binomial GLM to predict number of seedlings as a function of site
seedlings.mod <- glm.nb(tot.seedlings ~ Site, data = Just_seedlings)
summary(seedlings.mod)

## 
## Call:
## glm.nb(formula = tot.seedlings ~ Site, data = Just_seedlings, 
##     init.theta = 2.995739929, link = log)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.85486    0.19860  19.410  < 2e-16 ***
## SitePFCA 1   1.89803    0.27045   7.018 2.25e-12 ***
## SitePFCA 2   0.91920    0.27142   3.387 0.000707 ***
## SitePFCA 3   0.07303    0.27348   0.267 0.789434    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(2.9957) family taken to be 1)
## 
##     Null deviance: 111.465  on 38  degrees of freedom
## Residual deviance:  40.926  on 35  degrees of freedom
## AIC: 428.26
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  2.996 
##           Std. Err.:  0.672 
## 
##  2 x log-likelihood:  -418.256

Anova(seedlings.mod, type = "3")

## Analysis of Deviance Table (Type III tests)
## 
## Response: tot.seedlings
##      LR Chisq Df Pr(>Chisq)    
## Site   70.539  3  3.272e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Contrasts
standardize <- emmeans::emmeans(seedlings.mod, "Site")
emmeans::contrast(standardize, "trt.vs.ctrl")

##  contrast    estimate    SE  df z.ratio p.value
##  PFCA 1 - TP    1.898 0.270 Inf   7.018  <.0001
##  PFCA 2 - TP    0.919 0.271 Inf   3.387  0.0021
##  PFCA 3 - TP    0.073 0.273 Inf   0.267  0.9712
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: dunnettx method for 3 tests

## 95% confidence intervals 

emmeans::emmeans(seedlings.mod, "Site", type = "response")

##  Site   response    SE  df asymp.LCL asymp.UCL
##  TP         47.2  9.38 Inf      32.0      69.7
##  PFCA 1    315.1 57.84 Inf     219.9     451.5
##  PFCA 2    118.4 21.90 Inf      82.4     170.1
##  PFCA 3     50.8  9.55 Inf      35.1      73.4
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

## Make a new dataframe to hold the model estimates and 95% confidence interval bounds

response <- c( 315.1, 118.4, 50.8, 47.2)
lower <- c(219.9,  82.4,  35.1, 32.0) 
upper <- c(451.5, 170.1,  73.4,   69.7)
site <- c("PFCA 1", "PFCA 2", "PFCA 3", "TP")

seedlings.mod.est <- data.frame(site, response, lower, upper)


#Re-order the levels
Just_seedlings$Site <- factor(Just_seedlings$Site, levels=c("PFCA 1", "PFCA 2", "PFCA 3", "TP"))

### Plot number of seedlings vs. site

ggplot(Just_seedlings, aes(x= Site, y = tot.seedlings, color = Site))+
  geom_jitter(width =0.2, size = 2, show.legend = FALSE)+
  geom_point(data = Just_seedlings_mean, aes(x = Site, y = mean.seedlings),size = 3.5, color = "black")+
  geom_errorbar(data = seedlings.mod.est, aes(x= site, y = response, ymin = lower, ymax = upper), width = 0.15, cex = 1, color = "black")+
  theme_classic()+
  labs(x= "Site", y = "Number of Seedlings")+
  theme(text=element_text(size=18), legend.key.size=unit(1, "cm"), plot.title = element_text(hjust = 0.5))+
   scale_color_manual(name = "Site", breaks = c("PFCA 1", "PFCA 2", "PFCA 3", "TP"),
                         labels = c("Young", "Middle", "Old", "Remnant"),
                         values = c("#1f78b4", "#a6cee3", "#b2df8a","#33a02c"))+
   scale_x_discrete(name = "Site", breaks = c("PFCA 1", "PFCA 2", "PFCA 3", "TP"),
                         labels = c("Young", "Middle", "Old", "Remnant"))+  
  geom_vline(xintercept = 3.5, linetype = "longdash")+
  annotate("text", x = 1 , y = 850, label = "*", size = 9.5)+ 
  annotate("text", x = 2 , y = 300, label = "*", size = 9.5)+
  ylim(0, 1000)

Diversity

Summary stats

- Above - Total Species Richness

## # A tibble: 39 × 3
## # Groups:   Site [4]
##    Site   Transect Species_richness
##    <fct>  <fct>               <int>
##  1 TP     1                      34
##  2 TP     10                     37
##  3 TP     2                      30
##  4 TP     3                      31
##  5 TP     4                      27
##  6 TP     5                      30
##  7 TP     6                      29
##  8 TP     8                      38
##  9 TP     9                      43
## 10 PFCA 1 1                      28
## # ℹ 29 more rows

## 
## Call:
## glm(formula = Species_richness ~ Site, family = "poisson", data = Uni_Above_Diversity)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 3.503219   0.057831  60.576  < 2e-16 ***
## SitePFCA 1  0.121122   0.077532   1.562    0.118    
## SitePFCA 2  0.314493   0.074447   4.224  2.4e-05 ***
## SitePFCA 3  0.002338   0.079671   0.029    0.977    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 66.379  on 38  degrees of freedom
## Residual deviance: 40.952  on 35  degrees of freedom
## AIC: 261.23
## 
## Number of Fisher Scoring iterations: 4

## Analysis of Deviance Table (Type III tests)
## 
## Response: Species_richness
##      LR Chisq Df Pr(>Chisq)    
## Site   25.427  3  1.257e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##  contrast    estimate     SE  df z.ratio p.value
##  PFCA 1 - TP  0.12112 0.0775 Inf   1.562  0.2796
##  PFCA 2 - TP  0.31449 0.0744 Inf   4.224  0.0001
##  PFCA 3 - TP  0.00234 0.0797 Inf   0.029  0.9997
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: dunnettx method for 3 tests

##  Site   rate   SE  df asymp.LCL asymp.UCL
##  TP     33.2 1.92 Inf      29.7      37.2
##  PFCA 1 37.5 1.94 Inf      33.9      41.5
##  PFCA 2 45.5 2.13 Inf      41.5      49.9
##  PFCA 3 33.3 1.82 Inf      29.9      37.1
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

- Above - Total Species Richness

## # A tibble: 39 × 3
## # Groups:   Site [4]
##    Site   Transect Species_richness
##    <fct>  <fct>               <int>
##  1 TP     1                      34
##  2 TP     10                     37
##  3 TP     2                      30
##  4 TP     3                      31
##  5 TP     4                      27
##  6 TP     5                      30
##  7 TP     6                      29
##  8 TP     8                      38
##  9 TP     9                      43
## 10 PFCA 1 1                      28
## # ℹ 29 more rows

## 
## Call:
## glm(formula = Species_richness_native ~ Site, family = "poisson", 
##     data = Uni_Above_Diversity_native)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  3.44113    0.05965  57.684  < 2e-16 ***
## SitePFCA 1  -0.09123    0.08407  -1.085  0.27786    
## SitePFCA 2   0.20953    0.07846   2.670  0.00757 ** 
## SitePFCA 3  -0.06014    0.08343  -0.721  0.47102    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 58.191  on 38  degrees of freedom
## Residual deviance: 39.489  on 35  degrees of freedom
## AIC: 253.49
## 
## Number of Fisher Scoring iterations: 4

## Analysis of Deviance Table (Type III tests)
## 
## Response: Species_richness_native
##      LR Chisq Df Pr(>Chisq)    
## Site   18.703  3   0.000315 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##  contrast    estimate     SE  df z.ratio p.value
##  PFCA 1 - TP  -0.0912 0.0841 Inf  -1.085  0.5567
##  PFCA 2 - TP   0.2095 0.0785 Inf   2.670  0.0213
##  PFCA 3 - TP  -0.0601 0.0834 Inf  -0.721  0.7848
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: dunnettx method for 3 tests

##  Site   rate   SE  df asymp.LCL asymp.UCL
##  TP     31.2 1.86 Inf      27.8      35.1
##  PFCA 1 28.5 1.69 Inf      25.4      32.0
##  PFCA 2 38.5 1.96 Inf      34.8      42.5
##  PFCA 3 29.4 1.71 Inf      26.2      33.0
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

- SB - Total Species Richness

## # A tibble: 547 × 4
## # Groups:   Site, Transect [46]
##    Site    Transect SPP6   Number_Seedlings
##    <chr>   <fct>    <chr>             <dbl>
##  1 CONTROL 1        DIGISC                6
##  2 CONTROL 1        SETFAB                4
##  3 CONTROL 10       ALOCAR                1
##  4 CONTROL 10       CONCAN                1
##  5 CONTROL 3        ALOCAR                1
##  6 CONTROL 6        CONCAN                1
##  7 CONTROL 6        DIGISC                1
##  8 CONTROL 7        OXADIL                1
##  9 CONTROL 7        PANDIC                1
## 10 CONTROL 8        TAROFF                1
## # ℹ 537 more rows

## # A tibble: 40 × 3
## # Groups:   Site [4]
##    Site  Transect Species_richness
##    <fct> <fct>               <int>
##  1 TP    1                      14
##  2 TP    10                     13
##  3 TP    2                      11
##  4 TP    3                      11
##  5 TP    4                      12
##  6 TP    5                      10
##  7 TP    6                       4
##  8 TP    7                      12
##  9 TP    8                       9
## 10 TP    9                      15
## # ℹ 30 more rows

## 
## Call:
## glm(formula = Species_richness ~ Site, family = "poisson", data = Uni_Seedling_Diversity)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.40695    0.09492  25.359   <2e-16 ***
## SitePFCA 1   0.22494    0.12729   1.767   0.0772 .  
## SitePFCA 2   0.30775    0.12503   2.461   0.0138 *  
## SitePFCA 3   0.20312    0.12791   1.588   0.1123    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 53.060  on 39  degrees of freedom
## Residual deviance: 46.581  on 36  degrees of freedom
## AIC: 230.36
## 
## Number of Fisher Scoring iterations: 4

## Analysis of Deviance Table (Type III tests)
## 
## Response: Species_richness
##      LR Chisq Df Pr(>Chisq)  
## Site   6.4783  3    0.09052 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##  contrast    estimate    SE  df z.ratio p.value
##  PFCA 1 - TP    0.225 0.127 Inf   1.767  0.1919
##  PFCA 2 - TP    0.308 0.125 Inf   2.461  0.0382
##  PFCA 3 - TP    0.203 0.128 Inf   1.588  0.2674
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: dunnettx method for 3 tests

##  Site   rate   SE  df asymp.LCL asymp.UCL
##  TP     11.1 1.05 Inf      9.22      13.4
##  PFCA 1 13.9 1.18 Inf     11.77      16.4
##  PFCA 2 15.1 1.23 Inf     12.87      17.7
##  PFCA 3 13.6 1.17 Inf     11.50      16.1
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

- SB - Native Species Richness

## # A tibble: 40 × 3
## # Groups:   Site [4]
##    Site  Transect Species_richness_native
##    <fct> <fct>                      <int>
##  1 TP    1                             12
##  2 TP    10                            11
##  3 TP    2                             10
##  4 TP    3                             10
##  5 TP    4                             11
##  6 TP    5                              9
##  7 TP    6                              4
##  8 TP    7                             11
##  9 TP    8                              8
## 10 TP    9                             15
## # ℹ 30 more rows

## 
## Call:
## glm(formula = Species_richness_native ~ Site, family = "poisson", 
##     data = Uni_SB_Diversity_native)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.313e+00  9.950e-02  23.241   <2e-16 ***
## SitePFCA 1  -3.015e-02  1.418e-01  -0.213    0.832    
## SitePFCA 2   8.947e-11  1.407e-01   0.000    1.000    
## SitePFCA 3   4.832e-02  1.391e-01   0.347    0.728    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 55.512  on 39  degrees of freedom
## Residual deviance: 55.189  on 36  degrees of freedom
## AIC: 227.08
## 
## Number of Fisher Scoring iterations: 4

## Analysis of Deviance Table (Type III tests)
## 
## Response: Species_richness_native
##      LR Chisq Df Pr(>Chisq)
## Site  0.32366  3     0.9555

##  contrast    estimate    SE  df z.ratio p.value
##  PFCA 1 - TP  -0.0302 0.142 Inf  -0.213  0.9820
##  PFCA 2 - TP   0.0000 0.141 Inf   0.000  1.0000
##  PFCA 3 - TP   0.0483 0.139 Inf   0.347  0.9504
## 
## Results are given on the log (not the response) scale. 
## P value adjustment: dunnettx method for 3 tests

##  Site   rate   SE  df asymp.LCL asymp.UCL
##  TP     10.1 1.00 Inf      8.31      12.3
##  PFCA 1  9.8 0.99 Inf      8.04      11.9
##  PFCA 2 10.1 1.00 Inf      8.31      12.3
##  PFCA 3 10.6 1.03 Inf      8.76      12.8
## 
## Confidence level used: 0.95 
## Intervals are back-transformed from the log scale

PCoA

## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = everything_mat_hell_bray ~ Site * Type, data = everything_mat_labs, method = "bray", p.adjust = "bonferroni")
##            Df SumOfSqs      R2       F Pr(>F)    
## Site        3   10.843 0.31373 26.5607  0.001 ***
## Type        2    5.650 0.16347 20.7590  0.001 ***
## Site:Type   6    3.781 0.10940  4.6309  0.001 ***
## Residual  105   14.288 0.41341                   
## Total     116   34.561 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

## $parent_call
## [1] "everything_mat_hell_bray ~ Site * Type , strata = Null , permutations 999"
## 
## $`PFCA 1_vs_PFCA 2`
##           Df SumOfSqs      R2       F Pr(>F)    
## Site       1   1.9875 0.15857 16.8704  0.001 ***
## Type       2   3.4405 0.27449 14.6020  0.001 ***
## Site:Type  2   0.7445 0.05939  3.1596  0.001 ***
## Residual  54   6.3616 0.50755                   
## Total     59  12.5340 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $`PFCA 1_vs_PFCA 3`
##           Df SumOfSqs      R2       F Pr(>F)    
## Site       1   5.1421 0.30178 36.3215  0.001 ***
## Type       2   2.8375 0.16653 10.0215  0.001 ***
## Site:Type  2   1.4148 0.08303  4.9968  0.001 ***
## Residual  54   7.6448 0.44866                   
## Total     59  17.0392 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $`PFCA 1_vs_TP`
##           Df SumOfSqs      R2       F Pr(>F)    
## Site       1   6.0682 0.37289 48.6998  0.001 ***
## Type       2   2.3205 0.14259  9.3114  0.001 ***
## Site:Type  2   1.5299 0.09401  6.1391  0.001 ***
## Residual  51   6.3548 0.39050                   
## Total     56  16.2733 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $`PFCA 2_vs_PFCA 3`
##           Df SumOfSqs      R2       F Pr(>F)    
## Site       1   2.6237 0.16259 17.8593  0.001 ***
## Type       2   4.4234 0.27412 15.0549  0.001 ***
## Site:Type  2   1.1567 0.07168  3.9367  0.001 ***
## Residual  54   7.9331 0.49161                   
## Total     59  16.1368 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $`PFCA 2_vs_TP`
##           Df SumOfSqs      R2       F Pr(>F)    
## Site       1   3.8134 0.24391 29.2764  0.001 ***
## Type       2   3.6334 0.23240 13.9473  0.001 ***
## Site:Type  2   1.5447 0.09880  5.9296  0.001 ***
## Residual  51   6.6430 0.42489                   
## Total     56  15.6345 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $`PFCA 3_vs_TP`
##           Df SumOfSqs      R2       F Pr(>F)    
## Site       1   2.1080 0.13796 13.5635  0.001 ***
## Type       2   4.0498 0.26505 13.0290  0.001 ***
## Site:Type  2   1.1957 0.07825  3.8468  0.001 ***
## Residual  51   7.9262 0.51874                   
## Total     56  15.2797 1.00000                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## attr(,"class")
## [1] "pwadstrata" "list"

Indicator Species Analysis

Trajectory Analysis

## 
## Call:
## lm(formula = Directionality ~ Site_labels, data = trajectory.stats2)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.10512 -0.02697 -0.00408  0.04412  0.07062 
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.470789   0.008716  54.016  < 2e-16 ***
## Site_labelsYoung  -0.034596   0.012014  -2.880  0.00476 ** 
## Site_labelsMiddle -0.016128   0.012014  -1.342  0.18214    
## Site_labelsOld    -0.010394   0.012014  -0.865  0.38877    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04529 on 113 degrees of freedom
## Multiple R-squared:  0.07287,    Adjusted R-squared:  0.04826 
## F-statistic: 2.961 on 3 and 113 DF,  p-value: 0.03531

## Anova Table (Type III tests)
## 
## Response: Directionality
##             Sum Sq  Df   F value  Pr(>F)    
## (Intercept) 5.9843   1 2917.7332 < 2e-16 ***
## Site_labels 0.0182   3    2.9606 0.03531 *  
## Residuals   0.2318 113                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

##  contrast         estimate    SE  df t.ratio p.value
##  Young - Remnant   -0.0346 0.012 113  -2.880  0.0135
##  Middle - Remnant  -0.0161 0.012 113  -1.342  0.4009
##  Old - Remnant     -0.0104 0.012 113  -0.865  0.6994
## 
## P value adjustment: dunnettx method for 3 tests

Ternary Plot

## # A tibble: 21,177 × 5
## # Groups:   Site, Transect [39]
##    Site   Transect Type        SPP6   Abundance
##    <chr>  <fct>    <chr>       <chr>      <dbl>
##  1 PFCA 1 1        Aboveground ABUTHE         0
##  2 PFCA 1 1        Seed_Bank   ABUTHE         0
##  3 PFCA 1 1        Seed_Rain   ABUTHE         0
##  4 PFCA 1 10       Aboveground ABUTHE         0
##  5 PFCA 1 10       Seed_Bank   ABUTHE         0
##  6 PFCA 1 10       Seed_Rain   ABUTHE         0
##  7 PFCA 1 2        Aboveground ABUTHE         0
##  8 PFCA 1 2        Seed_Bank   ABUTHE         0
##  9 PFCA 1 2        Seed_Rain   ABUTHE         0
## 10 PFCA 1 3        Aboveground ABUTHE         0
## # ℹ 21,167 more rows

—-

Multipaned figures

## Warning in as_grob.default(plot): Cannot convert object of class list into a
## grob.